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Project  Abstract  (Original  Proposal) 


The  focus  of  this  research  program  is  to  combine  a  recently  developed  capability  for 
hybrid  large-eddy  /  Reynolds-averaged  Navier-Stokes  (LES/RANS)  simulations  with  a  novel 
immersed-boundary  method  to  predict  the  effects  of  different  techniques  for  forcing  boundary- 
layer  transition  and  mitigating  flow  separation  due  to  shock  /  boundary  layer  interactions.  Hybrid 
LES/RANS  methods  can  naturally  account  for  time-dependent  phenomena,  such  as  unsteady 
shock  motion,  and  have  displayed  an  ability  to  predict  the  rapid  recovery  of  a  turbulent  boundary 
layer  following  separation  and  re-attachment.  Because  the  near-wall  boundary-layer  structure  is 
modeled  through  RANS  concepts,  no  restrictions  in  Reynolds  number  are  required.  The  use  of 
immersed-boundary  techniques  promises  significant  economy  in  the  number  of  mesh  points 
required  to  simulate  complex  objects  such  as  boundary-layer  trip  arrays  and  perforated-plate 
bleed  systems.  Furthermore,  they  can  be  coupled  with  structural  response  models  in  a 
straightforward  manner,  thus  enabling  simulation  of  passive  control  methods  based  on  aero- 
elastic  flap  deflections.  The  proposed  research  program  is  focused  in  three  areas:  (1.) 
development  of  immersed-boundary  techniques  for  compressible,  turbulent  flows  at  high 
Reynolds  number;  (2.)  development  of  strategies  for  coupling  structural  response  models  with 
immersed-boundary  movement;  (3.)  simulation  of  well-documented  experiments  involving 
boundary-layer  trip  effects,  active  bleed  control  of  shock  /  boundary  layer  interactions,  and 
passive  control  of  shock  /  boundary  interactions  using  ‘meso  flap’  arrays;  and  (4.)  use  of  data 
from  hybrid  LES/RANS  simulations  to  guide  development  of  engineering-level  models  for  the 
effects  of  boundary-layer  control  devices. 
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1.  Background  and  Objectives 

The  control  of  shock-induced  separation  within  a  supersonic  inlet  is  critical  in  achieving 
the  desired  total  pressure  recovery  and  uniformity  of  the  flow  entering  the  compressor.  Off- 
design  responses,  such  as  inlet  unstart,  are  often  the  result  of  the  growth  of  large  regions  of 
separated  flow,  generated  initially  through  shock  interactions.  Control  of  shock-induced 
separation  within  the  inlet  often  involves  the  use  of  passive  or  active  bleed  systems  and  slots  to 
skim  off  the  boundary  layer  [1-5].  Control  techniques  based  on  aero-elastic  deformation  of  thin 
flaps  (UIUC’s  ‘mesoflap’  concept  and  others)  have  been  proposed  [6-8],  Micro  vortex 
generators  (micro  VGs)  are  currently  being  considered  as  a  means  of  energizing  a  boundary  layer 
through  the  introduction  of  axial  vorticity.  [9,10]  The  proper  implementation  of  such  control 
systems  requires  an  understanding  of  their  effects  at  on  and  off-design  conditions. 
Computational  fluid  dynamics  can  potentially  play  an  important  role  in  the  design  and 
implementation  of  boundary-layer  control  techniques,  but  the  direct  rendering  of  large  numbers 
of  such  control  devices  within  the  framework  of  a  complete  engine  simulation  is  currently 
infeasible.  Such  effects  must  be  modeled  at  some  level,  and  in  the  absence  of  detailed 
experimental  data,  there  is  little  available  to  guide  the  modeling. 

The  design  of  active  or  passive  bleed  systems  for  high-speed  inlets  is  an  area  in  which 
experimentation  and  empirical  correlations  have  played  a  dominant  role. [3,4]  Most  efforts  in 
simulation  of  boundary  layer  bleed  for  separation  control  have  concentrated  on  single  slot  or  hole 
effects  [11-14],  With  the  exception  of  a  few  three-dimensional  studies  involving  resolution  of 
individual  holes  using  overset-grid  techniques  [13,14],  most  studies  have  assumed  two- 
dimensional  flow,  and  all  have  utilized  Reynolds-averaged  Navier-Stokes  (RANS)  models  or 
simpler  strategies.  Given  that  fact  that  shock  /  boundary  layer  interactions  tend  to  be  unsteady 
on  a  large  scale,  and  that  local  pressure  differences  can  lead  to  periodic  blowing  /  suction  even  in 
“active”  control  devices  [5],  it  appears  that  accounting  for  the  time-dependent  nature  of  the 
control  approach  and  its  interaction  with  the  flow  may  be  critical  to  achieving  better  predictions. 
The  geometric  complexity  of  bleed  systems  complicates  the  application  of  techniques  such  as 
large-eddy  simulation  (LES),  as  resolution  of  features  such  as  individual  holes  would  normally 
require  a  large  number  of  clustered  mesh  cells.  The  number  of  mesh  points  used  could  become 
excessively  large,  and  the  accuracy  of  the  LES  approach  could  diminish  due  to  the  presence  of 
irregular  mesh  cells. 

Passive  control  of  shock  /  boundary  layer  interactions  (SBLIs)  through  the  use  of  vortex 
generators  embedded  within  the  boundary  layer  [9,10]  is  being  considered  as  an  alternative  to 
traditional  bleed  systems.  In  these  concepts,  micro  vortex  generators  are  embedded  in  regular 
arrays  within  the  boundary  layer  upstream  of  a  shock-impingement  point.  The  devices  modify 
the  structure  of  the  turbulent  boundary  layer,  energizing  to  the  point  that  it  can  better  resist  the 
shock-induced  adverse  pressure  gradient.  The  effects  of  flow-aligned,  regular  arrays  of  such 
devices  may  be  simulated  by  gridding  a  few  objects  in  detail  and  by  imposing  periodic  boundary 
conditions.  However,  this  approach  neglects  the  effects  of  possible  three-dimensionality,  which 
may  result  from  sidewall  shock  /  boundary  layer  interactions.  The  accurate  prediction  of  the 
performance  of  micro  VGs  may  also  require  a  time-dependent  method  of  analysis,  due  to  the 
need  to  capture  the  dynamics  of  the  generated  vortical  structures  and  their  effects  on  the  shock  / 
boundary  layer  interaction. 
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The  purpose  of  this  study  is  to  develop  a  different  framework  for  conducting  detailed 
time-dependent  simulations  of  the  effects  of  control  devices,  such  as  bleed  through  perforated 
plates,  aerolastically-deformable  ‘mesoflaps’,  and  micro  vortex  generators,  on  shock  wave  / 
boundary  layer  interactions.  The  approach  will  combine  hybrid  large-eddy  /  Reynolds-averaged 
(LES/RANS)  techniques  for  flowfield  modeling  [15,16]  with  a  novel  immersed-boundary  (IB) 
concept  [17]  for  geometric  rendering  of  control  devices.  The  use  of  LES/RANS  techniques  in 
this  scope  implies  that  the  effect  of  the  control  techniques  on  the  development  of  larger  turbulent 
eddies  should  be  captured  directly  for  improved  accuracy.  The  LES/RANS  concepts  developed 
at  NCSU  involve  the  use  of  flow-dependent  blending  functions  to  shift  the  closure  from  a  RANS 
description  near  solid  surfaces  to  a  LES  description  in  the  outer  portion  of  the  boundary  layer. 
Coherent  structure  growth  is  sustained  through  the  application  of  recycling  /  re-scaling  strategies 
[15].  As  the  near-wall  behavior  is  modeled  using  RANS,  the  total  number  of  mesh  points  needed 
is  significantly  less  than  for  a  wall-resolved  LES  at  the  same  Reynolds  number. 

The  use  of  immersed-boundary  (IB)  techniques  [18-20]  in  this  scope  promises  significant 
economy  in  terms  of  the  overall  number  of  grid  points  required,  as  fine  details  of  the  flow  near 
(or  within)  control  devices  are  modeled,  rather  than  directly  simulated.  The  hope  is  that  the 
geometric  influences  of  the  control  technique  on  the  surrounding  flow  are  captured  properly. 
Our  prior  work  [17]  for  incompressible  flows  has  shown  that  this  can  be  the  case  in  general,  even 
at  high  Reynolds  number.  Another  advantage  of  an  IB  method,  as  opposed  to  direct  rendering,  is 
that  feature  deformation  under  time-dependent  aero-elastic  loads  can  be  captured  simply  and 
directly  without  recourse  to  grid  re-generation  or  adaptation.  A  further  advantage  is  that  changes 
in  the  design  or  placement  or  number  of  control  devices  can  be  accomplished  in  a  trivial  manner, 
by  simply  re-generating  the  underlying  CAD  surface  descriptions  of  the  objects.  The  IB 
procedures  are  not  dependent  on  whether  the  turbulence  closure  is  RANS  or  hybrid  LES/RANS 
and  can  be  incorporated  into  production  level  RANS  codes  as  a  means  of  directly  simulating  the 
large-scale  effects  of  particular  control  methods. 

The  remainder  of  this  report  outlines  the  numerical  methods  used  in  Section  2  and  the 
turbulence  closure  methods  employed  in  Section  3.  A  complete  description  of  an  IB  method 
valid  for  compressible,  turbulent  flows  is  provided  in  Section  4,  and  the  extension  of  the 
methodology  to  aero-elastic  deformation  of  surfaces  as  described  by  beam-bending  theory  is 
discussed  in  Section  5.  The  results  of  the  investigation  are  presented  in  Section  6.  Three  major 
classes  of  control  devices  for  shock  /  boundary  layer  interactions  are  considered:  micro  VGs 
(Section  6.2),  bleed  arrays  (Section  6.3),  and  mesoflaps  (Section  6.4).  In  addition,  results  from 
other  studies  that  employ  the  developed  IB  method  to  simulate  laminar  flows  over  boundary 
layer  trip  devices  (Section  6.5)  and  shock-generator  effects  in  a  3-D  shock  /  boundary  layer 
interaction  (Section  6.6)  are  described  in  brief.  Conclusions  of  the  study  are  presented  in 
Section  7,  and  directions  for  future  work  are  discussed  in  Section  8. 

2.  Numerical  Methods 

Simulations  were  performed  using  NCSU’s  REACTMB  code.  In  REACTMB,  the 
governing  equations  are  discretized  in  a  finite -volume  manner.  Inviscid  fluxes  are  discretized 
using  Edwards’  LDFSS  scheme  [21],  while  viscous  and  diffusive  fluxes  are  discretized  using 
second-order  central  differences.  The  LDFSS  scheme  is  extended  to  second  or  higher-order 
spatial  accuracy  using  the  piecewise  parabolic  method  (PPM)  [22],  as  described  later.  A  dual¬ 
time  stepping  implicit  method  is  used  to  advance  the  equations  in  time.  At  each  time  step,  a 
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Crank-Nicholson  discretization  of  the  equations  is  solved  to  a  prescribed  tolerance  using  a  sub¬ 
iteration  procedure.  The  matrix  system  resulting  from  the  linearization  of  the  equation  system  is 
approximately  solved  using  a  planar  relaxation  procedure  at  each  sub-iteration.  To  enhance 
computational  efficiency,  matrix  elements  are  held  fixed  over  the  duration  of  the  sub-iterations. 

The  PPM  reduces  to  a  fourth-order  central  differencing  scheme  for  sufficiently  smooth 
data,  enhancing  its  ability  to  capture  small-scale  features.  A  cell-by-cell  limiting  procedure 
reduces  the  order  of  accuracy  in  the  vicinity  of  local  extrema  in  the  solution.  A  general  interface 

flux  formula  may  be  written  as  Fi+U2(VLi+U2,VRM/2) ,  with  left-and-right  state  interpolations 
VLi+ 1/2  and  VRM/2  determined  by  the  step-by-step  procedure  described  below. 


Step  1:  initially  set  left-and-right  state  values  to  fourth-order  averaging  operators.  If 
unmodified,  these  provide  a  fourth-order  central  difference  approximation  to  the  flux 


difference  Fi+U2  -  Fi_U2 : 


i+l/2  —  ^R,i+U2  ~  2  <y,+yM)~<yH2+yM) 


in 


Step  2:  conduct  cell-by-cell  limiting  to  enforce  monotonicity: 

if  sgn[(Zu+„2- W  -iV„2)]  =  -l,  then 


^Z,  t+1/2  —  1-1/2  —  ^ 


else 


c  =  V  -V 

y  L,  1+1/2  y  R,  1-1/2 

D  =  6. 0K-i(F£,+1/2+Ffi,-1/2)] 
if  (DC  >  CC)  then 

V  =  3V -2V 

y  R,i- 1/2  Jy  i  Z, 1+1/2 

elseif  (-CC  >  DC)  then 

V  =3V -2V 

y  L, ;+l/2  Jyi  ^y  R,i- 1/2 

end  if 
end  if 


(2) 


The  first  ‘if  block  resets  the  interpolation  function  to  a  constant  if  Vi  is  a  local  maximum  or 

minimum.  The  second  ‘if  block  resets  either  the  left-state  value  at  interface  i  + 1/2  or  the  right- 
state  value  at  interface  i  - 1  /  2  so  that  the  interpolation  parabola  that  connects  the  interface  states 
with  the  state  at  the  cell  center  is  monotonically  increasing  or  decreasing.  It  is  also  possible  to 
enforce  physical  constraints  (such  as  positive  temperatures,  densities,  and  mass  fractions)  at  the 
cell  interfaces  by  a  similar  cell-by-cell  resetting  algorithm.  The  PPM  requires  a  seven  point 
stencil  in  each  coordinate  direction,  and  the  reconstruction  procedures  are  applied  to  the 

primitive-variable  vector  V  =  [p, u, v,  w, T, k, co]T . 


7 


3.  Turbulence  Closure 


3.1.  Menter’s  SST  Model 


The  RANS  turbulence  model  used  in  this  investigation  is  Menter’s  hybrid 
k-col  k-s shear-stress  transport  (SST)  model. [23]  The  transport  equations  for  the  turbulence 
kinetic  energy,  k ,  and  the  specific  dissipation  rate,  co  for  Menter’s  model  are  given  by 

d{pk)  d(pkuj )  2  a.  ,  d 

—  l  +  — - —  =  u.S  -  B  pcok-v - 

dt  dx:  dx : 


ox, 


(3) 


dt 


dx. 


dxj 


dco 

dxj 


1  dk  dco  ... 
+  2(1-/^  )paai2  — — — (4) 


cfcc .  dxj 


where  S  is  taken  to  be  the  magnitude  of  the  vorticity  vector  and  <j)  represents  any  constant  in 
Menter’s  model  [crk  averaged  according  to  the  blending  function  Fx 

<t>  =  FA  +  (l  -  A,  V2  (5) 

The  <j)x  constants  are  from  the  k  -  co  model: 

<rki  =  0-85,  <7ml  =  0.5  ,  px  =  0.075 


jg  =  0.09 ,  k  =  0.41 ,  yx  =  - 


A 


F  FF 


and  the  constants  are  from  the  standard  k-  s  model: 


< 7k2  =1.0,  Gal  =0.856,  P2  =0.0828 


(6) 

(7) 


p*  =0.09,  k  =  0.41 ,  y2  = 


A 


A’  V7 


500tA 


0.09at/  ’  d2co 


The  blending  function  Fx  is  defined  as 

Fx  =  tanh(arg^ )  ,  arg,  =  max 

V - -  "  ”  J 

where  d  is  the  distance  to  the  nearest  wall.  The  eddy  viscosity  is  defined  as: 

a,k 

v  — - - - 

‘  m&x{axo),  Ass/Q.F2  ) 

where  a,  =0.31,  Q  is  the  magnitude  of  the  vorticity  vector,  and  F2  is  another  blending  function 
given  by 


(8) 


(9) 


F2  =  tanh(arg2),  arg2  =  max 


4k  500t> 


v  0.09 cod  ’  d2co  j 


(10) 


This  ‘shear  stress  transport’  (SST)  modification  alters  the  turbulence  frequency  as  used  in  the 
eddy  viscosity  definition  in  the  vicinity  of  high  strain  rates.  The  general  effect  is  to  lower  the 
eddy  viscosity  and  thus  to  promote  the  growth  of  regions  of  separated  flow.  For  the  relatively 
weak  interactions  considered  in  this  study,  the  use  of  the  SST  modification  typically  results  in 
better  RANS  predictions.  The  Menter  BSL  model  is  used  as  the  RANS  component  of  the  hybrid 
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LES/RANS  model,  as  previous  studies  [15]  have  found  that  this  approach  provides  generally 
better  results  than  the  use  of  the  SST  model. 


3.2.  Hybrid  LES/RANS  Model 


The  LES/RANS  model  used  in  this  work  is  designed  to  transition  from  RANS  near  solid 
surfaces  to  LES  as  the  boundary  layer  shifts  from  its  logarithmic  form  to  its  wake-like  form.  The 
eddy  viscosity  as  used  in  the  turbulence  transport  equations  and  the  main  flow  equations  is 
defined  as  a  weighted  sum  of  a  RANS  description  and  a  mixed-scale  model  [24]: 

k  ^ 


m,  =  pvt  =  p\  r  —  +  (i  -  f> t,sGs 

CO  J 


(11) 


with 


v. 


=  CMSV2(q2)l,4A3,2,CM=0.06 


't,SGS~^M^  VV  J  i  '-'M  v'-v'v’  (12) 

An  estimate  of  the  subgrid  kinetic  energy  is  obtained  by  test-filtering  the  resolved-scale  velocity 
data: 

1 


4  =-(«* 


uky 


(13) 


The  blending  function  T  is  designed  to  shift  the  closure  from  a  RANS  description  near 
solid  surfaces  to  an  LES  description  in  the  outer  parts  of  the  boundary  layer  and  in  regions  of 
flow  separation.  The  blending  function  is  based  on  the  ratio  of  the  wall  distance  d  to  a  modeled 
form  of  the  Taylor  micro-scale: 


r  = 


l 


K 


1  -  tanh[  5(  . —  77  -1)-^] 


Vc, 


77 


(14) 


with  the  Taylor  micro-scale  defined  as 


A  =  Jv  ~j~C f 


CO 


(15) 


The  primary  advantage  of  this  form  is  that  the  location  of  the  RANS/LES  juncture  (defined  as 
T=0.5)  can  be  correlated  as  a  function  of  the  wall  coordinate  d+  =  urd  /  v  .  Substituting  the  log- 


law  expression  for  the  turbulence  frequency  co  =  uT  /(^[fA/cd)  into  Eq.  (15)  and  placing  the  result 
in  Eq.  (14),  it  can  be  shown  that  the  r=0.5  position  should  occur  at  d+  =  a 2  if  <j>=  0.0.  The 
constant  cq  then  can  be  used  control  the  position  of  the  juncture,  assuming  that  the  log-law 
scaling  is  maintained.  The  value  of  “5”  in  Eq.  (14)  controls  the  sharpness  of  the  transition,  and 
the  constant  (j)  can  be  used  to  shift  the  mid-point  of  the  blending  function.  In  this  work,  (j)  is  set 


to  tanh  ‘(0.98)  ,  so  that  the  balancing  position  (where 


=  1 )  corresponds  to  T  =0.99 


instead  of  T=0.5. 

To  determine  the  constant  ax  for  a  particular  inflow  boundary  layer,  the  following 
procedure  is  used.  First,  a  prediction  of  the  equilibrium  boundary  layer  is  obtained,  given  free- 
stream  properties,  a  specified  wall  condition  (adiabatic  or  isothermal)  and  a  value  for  the 
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boundary  layer  thickness,  from  Coles’  Law  of  the  Wall  /  Wake  along  with  the  Van  Driest 
transformation: 


u,„,  1  ,  ,  „  fl  .  ltn  d v  ,4.  urd 


^  =  -ln«)  +  C  +  2-shV(--),  d+w=- 
K  K  2  O  V 


with 


v4  = 


-i 


sin 


2  A2ul  -  B 

Jb^+aa1 


+  sin 


5 


V#2  +  4vf2 


(16) 


(17) 


^  Pr  M2  ,  B  = 

i+P  i/2(r  !)M2 

<  2  '  Tw 

2 

T 

1 

r„ 


An  initial  estimate  for  the  outer  extent  of  the  log  layer  is  defined  by  finding  the  value  of  d such 


that 


-In  OO  +  C 


A  ^  1 

vd  I  =  0.98 .  The  value  of  d+  =  urd  Iv  that  corresponds  to  this  value  of 


\UT  J 


dl  is  then  found  through  the  use  of  Walz’s  formula  for  the  static  temperature  distribution  within 
the  boundary  layer  [25]: 


+  -O'"1)  „2 


-M: 


f u * 


VWoo  J 


(18) 


Too  L  Too  Woo  2 

Since  the  kinematic  viscosity  v  is  a  function  of  temperature,  the  target  value  for  d+  will  differ 


significantly  from  d+  for  high  Mach  number  flows.  The  selection  of  the  constant  is  done  in  a 

pre-processing  step,  based  on  an  inputted  distribution  of  the  boundary  layer  thickness  versus 
streamwise  distance  obtained  from  a  RANS  solver.  Table  1  shows  streamwise  distributions  in 
the  model  constant  for  the  cases  described  later.  A  reasonable  estimate  of  the  value  of  model 
constant  can  be  obtained  by  considering  only  the  first  term  in  the  fit,  as  the  streamwise  interval  is 
generally  small  (~10  boundary  layer  thicknesses) 


Table  1:  Model  constant  versus  streamwise  distance 


case 

Model  constant  a, 

Babinsky,  et  al.  [10] 

a,(x)  =  27.35  +  39.02  x  - 19.21  x2 

Willis,  et  al.  [26] 

«j(x)  =  53.90  + 1 1.04  x 

Gefroh,  et  al.  [6] 

«j(x)  =  38.46  + 49.78  x 

A  key  to  the  blending  function  response  is  the  relative  insensitivity  of  the  turbulence 
frequency  to  changes  in  the  turbulence  kinetic  energy  and  the  eddy  viscosity.  This  means  that 
the  average  location  of  the  blending  function  (in  terms  of  its  distance  from  the  wall)  will  not  vary 
significantly.  The  dependence  of  the  turbulence  frequency  on  the  vorticity  S  means  that  the 
instantaneous  position  of  the  interface  will  vary  in  response  to  resolved-eddy  dynamics.  The 
calibration  procedure  described  above  is  specific  to  a  particular  inflow  boundary  layer  and  is  not 
guaranteed  to  adjust  properly  as  the  flow  transitions  from  one  equilibrium  state  to  another.  An 
improved  method  that  does  not  require  a  problem-specific  calibration  has  been  developed  and  is 
undergoing  evaluation  [27]  but  has  only  been  used  for  the  case  presented  in  Section  6.6. 
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3.3.  Recycling-Re-scaling  Method 


A  recycling  /  rescaling  technique,  applied  to  the  fluctuating  fields,  is  used  to  initiate  and 
sustain  turbulent  structures.  In  this  procedure,  fluctuations  in  the  fluid  properties  at  a  ‘recycle 
plane’  are  extracted  by  subtracting  the  instantaneous  profde  from  a  time-  and  span-averaged 
profile.  The  fluctuation  fields  are  then  rescaled  according  to  boundary  layer  similarity  laws  and 
superimposed  onto  a  RANS  mean  inflow  profde.  The  general  procedure  is  described  in  [15].  One 
modification  used  in  this  work  to  prevent  excessive  turbulence  energy  accumulation  in  the  outer 
part  of  the  boundary  layer  is  to  multiply  the  recycled  fluctuations  by  a  Klebanoff-type 
intermittency  function  [28]: 

inflow  =  ?rans  +  Fkieb  Cyc  J  C1  =  density,  velocity,  temperature 


inflow  =  (1  -  ^kieb  )</rans  +  ^kieb  4,'ecyc ;  q  =  turbulence  variables  ( 1 9) 

Fkkh=(\+(d/caebs)6yl 

with  Ckleb=  1.10.  This  procedure  also  ensures  that  the  RANS  ffee-stream  inflow  properties  are 


not  altered  outside  the  boundary  layer.  We  also  require  that  recycled  temperature  fluctuations  be 
no  greater  (in  magnitude)  than  allowed  by  Morkovin’s  hypothesis  of  negligible  total- temperature 
fluctuations.  From  Morkovin’s  hypothesis, 


CpT mork  —  [^kleb(M 


1  T7 

2  1  kleb 


RANS  U  recyc  +  VRANSVrecyc  +  WRANSVViecyc  )  + 


+  v'2  +  w'2 

recyc  recyc  recyc 


)] 


(20) 


and  the  recycled  temperature  fluctuation  is  limited  as  follows: 


Cycllm=^klebCycmin(1.0, 


kleb  recyc 


T 

mork 


F  T' 

kleb  recyc 


), 


^kleb^recyc  ^  ^ 


(21) 


The  density  fluctuation  is  determined  as  follows.  First,  a  provisional  value  of  the  pressure 
fluctuation  is  determined  from  the  recycled  density  fluctuation: 


prov 


■^[/^RANS^recyc  bm  ^kleb/^recyc(^RANS  ^recyc  ,j„  )1 


(22) 


The  pressure  fluctuation  is  then  limited  to  be  a  specified  multiple  of  the  pressure  (2%  in  this 
case): 

/Vov|lim  =  sgn(/?prov  )  min(|//rov  |,0.02/>RANS )  (23) 

and  a  corrected  value  of  the  density  fluctuation  is  determined  from  the  limited  pressure 
fluctuation: 


F recyc  (P  Pr°v  |  |jm  /^RANS^ ^recyc [.  )  KARANS  +  ^recyc|i;m  ) 


lim 


recyc  I  lim  } 


(24) 


A  further  modification  is  motivated  by  the  fact  that  quasi-periodic  inflow-generation 
methods  tend  to  ‘fix’  the  positions  of  longitudinal  vortical  structures  within  the  boundary  layer 
over  long  times.  This  problem  is  exacerbated  by  the  RANS  component  of  the  closure,  which 
attenuates  smaller-scale  near-wall  structures  that  might  otherwise  act  to  break  up  the  quasi- 
periodic  behavior.  One  way  to  avoid  this  is  by  shifting  the  plane  containing  the  fluctuation  fields 
in  the  spanwise  (Z)  direction  by  some  fraction  of  the  incoming  boundary  layer  thickness  8 .  This 
shift  is  determined  by  a  profile  of  shifting  distance  versus  time  that  is  produced  by  sampling 
Gaussian  distributions  corresponding  to  an  average  near-wall  streak  length  (3-5  8 )  and  the 
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shifting-distance  itself.  To  produce  this  profile,  the  average  near-wall  streak  length  (usually  3-5 
boundary  layer  thicknesses)  is  converted  to  an  increment  in  time  by  dividing  the  value  by  the 
convective  velocity  (taken  as  80%  of  freestream  velocity).  Standard  deviations  for  the  streak- 
length  and  spanwise-  shifting  Gaussian 
distributions  are  taken  as  1.0  8  and  as 
0.5  8  ,  respectively.  The  procedure 
results  in  a  random  sequence  of  shifts  at 
discrete  time  intervals.  This  discrete 
distribution  is  converted  to  a  continuous 
one  by  linearly  interpolating  between 
each  discrete  point.  A  typical 
distribution  of  shifting  distance 
(normalized  by  8  )  versus  normalized 
time  is  shown  in  Figure  1.  Ref.  [29] 
shows  that  this  procedure  does  act  to 
homogenize  the  time-averaged  flow  in 
the  spanwise  direction  while  not 
significantly  affecting  turbulent 
statistics.  The  ‘Z-shifting’  technique 
was  only  used  for  the  simulations  of  the 
Willis,  et  al.  [26]  and  Gefroh,  et  al.  [6] 
experiments. 


Figure  1:  Spanwise  shifting  distance  versus  time 


4.  Immersed  Boundary  Method 


In  the  present  work,  the  IB  method  developed  in  Choi  et  al.  [17]  is  extended  to  handle 
compressible,  turbulent  flows.  The  immersed  surface  is  generated  as  a  cloud  of  points  which  can 
be  structured  or  unstructured.  The  entire  flow  domain  is  then  classified  into  three  categories  of 
cells.  Cells  sufficiently  removed  from  the  immersed  boundary  are  termed  as  ‘field’  cells,  cells 
very  near  but  not  inside  the  immersed  object  are  ‘band’  cells,  and  cells  inside  the  immersed  body 
are  ‘interior’  cells.  To  perform  this  classification,  a  distance  function  to  the  nearest  surface  point 
for  all  cells  (within  a  ‘bounding  box’  of  the  IB)  is  computed  using  an  approximate  nearest 
neighbor-search  algorithm.  Concepts  from  computational  geometry  are  then  used  to  impose  an 
unambiguous  sign  to  the  distance  function  (thus  deciding  ‘inside’  or  ‘outside’  for  closed 
surfaces).  A  ‘direct  forcing’  approach  is  used  to  enforce  the  boundary  conditions  at  the  interior 
and  band  cells.  This  results  in  the  residual  form  of  the  equation  system  shown  below  which  is 
then  solved  implicitly,  coupled  with  exterior  cells,  by  use  of  sub-iteration  techniques. 


Rn+U  =  (q_  G(<p«+‘))£»+M  +  G(0"+1) 


n+ 1 ' 


V, 


n+l,l 


■V, 


n+l,l ' 


B,i 


At 


=  0 


(25) 


This  equation  represents  the  blending  of  the  Navier-Stokes  residual  with  a  source  term  that 
relaxes  the  primitive  variable  vector!7  =  [p,u,v,w,T,k,co]T  to  its  band-cell  values.  The  quantity 
G(0 )  is  a  sharp  Heaviside  function  (set  to  1  for  ‘band’  and  ‘interior’  cells  and  zero  otherwise), 
O  is  the  signed  distance  function,  and  /  is  a  sub-iteration  index.  Field  cells  are  defined  as  those 
with  G  =  0  and  0  >  0,  while  band  cells  are  those  with  G  =  1  and  0  <  0,  and  interior  cells  are 
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those  with  G  =  1  and  O  <  0.  The  Navier-Stokes  equations  are  solved  in  the  field  cells,  fluid 
properties  are  generally  held  fixed  in  the  interior  cells,  and  analytic  forms  for  the  fluid  properties 
(discussed  later)  are  prescribed  in  the  band  cells. 

4.1.  Cell  Classification  Procedure 

The  first  step  in  the  IB  method  is  to  define  an  immersed  body  as  a  collection  of  surface 
points.  This  can  be  done  using  a  CAD  format  or  through  other  means,  but  the  key  is  that  for 
each  separate  component  (separate  in  the  sense  that  different  components  may  move  at  different 
rates),  a  list  of  surface  points  and  outward-pointing  unit  normals  is  created.  The  next  step  is  to 
calculate  the  signed  distance  from  each  field  point  xk  to  the  nearest  surface  point  on  each  surface 

xsMk) .  This  is  accomplished  first  through  the  use  of  approximate  nearest-neighbor  searching 

techniques  [30],  which  return  the  unsigned  distance.  In  practice,  this  is  done  only  for  the 
number  of  field  points  that  are  within  a  “bounding  box”  surrounding  the  particular  surface,  as  it 
is  these  points  that  are  likely  to  be  influenced  immediately  by  the  body.  Distances  outside  the 
“bounding  box”  are  assigned  to  be  a  very  large  positive  number.  For  the  cases  presented  in  this 
report,  the  signed  distance  function  is  obtained  by  multiplying  the  unsigned  distance  with  the 
sign  of  the  dot  product  of  the  distance  vector  with  the  outward  normal: 

=  s§n((4  -*,,/(*) )•«.,/«  )x  I  ^  “*,,/(*)  I  (26) 

Other  approaches  [17]  can  be  used  for  more  complicated  CAD-based  objects. 

To  define  a  global  signed  distance  function  at  any  given  mesh  point  (0(4,0),  a  simple 

priority  rule  is  exercised.  First,  the  global  distance  function  is  initialized  to  a  large  number. 
Then,  the  global  signed  distance  function  at  a  particular  point  is  taken  as  the  minimum  of  the 
individual  signed  distance  functions  at  that  point: 

0(4,0  =  min,  (<I>,  (4,0)  (27) 

The  collections  of  points  that  comprise  the  surfaces  are  allowed  to  move  according  to  prescribed 
rate  laws.  Once  new  surface  positions  are  obtained,  the  preceding  and  following  steps  for 
developing  the  signed-distance  and  Heaviside  functions  have  to  be  repeated. 

The  Heaviside  function  G(O(x,0)  is  defined  to  be  one  for  points  just  outside  the 
immersed  body  and  within  the  immersed  body  and  is  zero  otherwise.  The  calculation  of  the 
Heaviside  function  is  initiated  by  first  initializing  G(O(4,0)  =  0  f°r  all  points xk .  Then,  given  a 

point 4  ,  if  <t> (4,0  >  0  and  if  any  0(4  ,t)  <  0,  where  41s  a  nearest-neighbor  of  xk ,  then 
<7(0(4, 0)  is  set  to  1.  If  0(4,0  <  0,  then  G(O(4,0)  Is  also  set  t°  1-  The  set  of  “nearest- 
neighbors”,  for  a  structured  grid,  is  generally  defined  as  the  26  cells  that  are  immediately 
adjacent  to  a  particular  mesh  cell  (i,j,k),  though  smaller  subsets  can  be  used.  As  mentioned 
before,  this  process  subdivides  all  of  the  mesh  cells  into 

a. )  field  points  (those  with  0(4, t)>  0  and  G(O(4,0)=0) 

b. )  band  points  (those  with  0(4, t)>  0  and  G(O(4,0)=0,  and  (28) 

c. )  interior  points  (those  with  0(4  >t)<  0  and  G(O(4,0)=l) 
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4.2.  Band  Cell  Interpolations 


To  extend  the  method  of  [17]  to  compressible  flows,  the  following  first-order  accurate 
closures  are  considered  for  the  fluid  properties  in  the  band  cells,  where  the  subscript  T  indicates 
properties  obtained  at  an  interpolation  point  located  along  the  normal  line  extending  outward 
from  the  nearest  surface  point  corresponding  to  the  band  cell  in  question  (discussed  later),  and 
the  subscript  “B”  indicates  the  band-cell. 

Ps  =  P(di ) 


f  j  A 


*BA 


’  US,i  ~  UT,i  idI ) 


ydI  j 


+  uNi(dj)c(p,dj,dB), 


(29) 


UnM,)  =  (uj(di)~usj)njni’  uTi(dj)  =  (iil(dl)-aSi)-aNi(dl) 

In  these  expressions,  n  is  the  normal  vector  at  the  closest  point  on  the  body  surface,  d  is  a 
distance  from  the  nearest  surface  point,  uS  j  is  the  velocity  at  the  nearest  surface  point,  and  A:  is  a 

power-law.  The  choice  of  k  allows  the  model  to  replicate  a  turbulent  velocity  profile  (£=1/7  or 
1/9)  or  a  laminar  profile  (k=  1).  To  obtain  the  temperature  distribution  near  the  surface,  Walz’s 
relation  for  the  temperature  distribution  within  a  compressible  boundary  layer  is  used  [25]: 


T(d,)  T(d,) 


■  + 


T  r(y-V) 

1 - — +  — - —Wnid,)] 

T(dj)  2yRT(d1)  Ta  , 


Y  J  Y 


V  dl  J 


Kr-i) 

2yRT(d,) 


[uT  j  {d  j )] 


(  d^ 

\d i  j 


2k 


or - 


=  1  +  - 


r(r~  1) 


K,(^)]2(1- 


(dg\ 

\di  J 


2k 


(30) 


)  (adiabatic  wall) 


T(dj)  2yRT(dj) 

In  this,  r  is  the  recovery  factor  and  [uT t(d ,  )]2  is  the  kinetic  energy  associated  with  the  tangential 


velocity  component  at  the  interpolation  point. 

The  function  c(p,dr,dB)  that  scales  the  normal  velocity  component  in  Eq.  (29)  is 
determined  by  enforcing  a  discrete  form  of  the  continuity  equation  at  each  band  cell.  In  the 
immersed  boundary  method  of  Choi  et  al.  [17],  we  adopt  a  locally-parallel  flow  assumption  so 
that  all  flow  properties  in  the  region  near  the  immersed  surface  are  functions  of  the  coordinate 
normal  to  the  surface  n.  Decomposing  the  velocity  into  tangential  components  uT  and  aT  (in 


directions Tx  and  T2)  and  normal  component un (in  direction  n),  we  write  the  steady  continuity 
equation  as: 


V  T  •  (^puT  )  ■+ 


8n 


=  0,  with  V  T  ■  (puT  )  = 


d{puT  )  d(puT  ) 


dTx 


-  +  - 


ar. 


(31) 


Applying  the  parallel-flow  assumption  and  introducing  the  functional  forms  / ( n )  and  g(n)  to 
describe  the  distributions  in  tangential  velocity  and  density  from  the  interpolation-point  location 
dj  to  the  wall  (ie.  uT  ( n )  =  uT  (d, ) / ( n ),  p{n)  =  p(dr )  f  ( n ) ),  we  have 


V T  •  ( p(dj)uT  (d j))f(n)g(n )  I  =  0 

on 


(32) 
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Now  consider  the  two  control  volumes  Qj  and  Q2  shown  in  Figure  2.  Discretizing  Eq.  (32)  over 
each  of  the  control  volumes  gives  the  following  system: 


Figure  2:  Schematic  of  control  volumes  used 
to  define  normal  velocity  component 


V  T  •  (p(dj  )uT  (d ,))f(d  +  )g(d  +  )  4- 


P(dj )u N (dI )  p(d B^u N  (d B) 

d  j  d  g 


=  0,  Q  =  Qj 


V r  •  {p{d,  )uT  (d ,))f  (d  -  )g(d  -  )  =  Qj  Q  =  Q2 


(33) 


d  r 


with 


1 


d  —  (d  j  -\-dB), 


1 


(34) 


d  =  —  d 


Upon  eliminating  the  common  factor Vr  ■(p(dI)uT(d 7)) ,  Eq.  (33)  can  be  solved  to  obtain  an 
expression  for  the  normal  velocity  component  at  the  band  cell,  uN  (d B  ) : 


Uff(dB)  Pidl) 


un(^i)  p(d  B) 


^/(J-)g(J-) 

dI 


“P  f(d  -  )g(d-)  +  f(d')g(d  *  )(1  -  df) 

v  a i  J 


(35) 


The  closure  is  completed  by  substituting  power-law  and  Walz  /  Crocco-type  relations  for 
f(n)  andg(n) : 


fin)  = 


(  Y 

'  n  > 


\di  j 


»  gin)  = 


T{n) 


(from  Eqs.  (29)  and  (30)),  (36) 
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by  assuming  that  the  pressure  gradient  normal  to  the  surface  is  zero,  and  by  evaluating  the 
expression  at  the  band-cell  distance  dB .  The  result,  for  an  adiabatic  wall,  is  given  as 

dn 


c(p,dj,dB)  = 


d~p ~ 

1  d, 

dI  dj 


{  dK  1 

k 

(  d  d  } 

d  = 

,  d~  = 

i(i+^r) 

^2  dI  j 

k  UI  J 

i=i+ 


r(r- 1) 


P  2  rRTid,) 


K,K)]2(1- 


j 


2k 


X 


■k=1+ irwh i«TAd,)f(i-(d-n 
p  2/RT(dj) 

■k = ' + r,  k  )t  O' -  (</* 

p  2/RT(dj) 


(37) 


The  turbulence  variables  in  the  band  cells  are  defined  as 

kB  =  u2 /jc;  ,  /(JC^/<dB) :  d+  >  10.934 


k(dj) 


dj 


,  coB  =  60vw  /(0.075d2B)  : 


d+  <  10.934  or  k  =  \ 


(38) 


=  uTdb  lvw,  uT  =|  uT  i (d, )  |  /(ln(d+)/ at  +  5.1)  (iterative  solution) 

To  arrive  at  this  form,  we  assume  equivalence  between  the  result  provided  by  the  power-law 
profile  and  the  law  of  the  wall  within  the  band  cells.  Some  of  the  cases  presented  later  use 
another  variant  of  this  procedure  that  replaces  the  pressure  extrapolation  (the  first  line  in  Eq. 
(29))  with  the  finite-volume  solution  of  the  continuity  equation  within  the  band  cells.  In  the 
formulation  of  the  continuity  equation,  the  mass  flux  is  set  to  zero  at  a  cell  interface  if  the 
distance  function  cD  changes  sign  across  the  interface.  The  continuity  equation  in  the  band  cells 
evolves  according  to  this  constraint  and  to  the  specified  velocity  and  temperature  fields  described 
above.  This  procedure  ensures  exact  mass  conservation  but  can  sometimes  lead  to  pressure 
oscillations  within  the  band  cells. 


4.3.  Determination  of  Information  at  the  Interpolation  Point 

The  preceding  developments  hinge  on  the  determination  of  flow  properties  q(dj )  at  a 
certain  distance  d:  away  from  the  surface.  Given  a  point  within  the  band  xk  and  a  list  of  nearest 
neighbors  to  that  point  3c, ,  a  merit  function  w,  is  defined  as  follows:  (see  Figure  3) 

Wj  =  =  - -  if  (X[  -xk)-h>  0,  and  wt  =  0  otherwise  (39) 

V(l*/ -**  I)2  “((*/  ~xk)-n))2  +e 
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In  this,  {xl  -  xk )  ■  h  is  the  projection  of  the  distance  from  xk  to  xl  in  the  direction  of  the  outward 
normal,  and  |  x,  -  xk  |  is  the  magnitude  of  the  distance  vector  itself.  If  point  x,  is  located  directly 
along  the  outward  normal  line  corresponding  to  band  point  xk ,  and  if  (xf  -  xk )  •  h  is  positive, 


Figure  3:  Schematic  determination  of  the  distance  dk  between  the  interpolation  point 
a  and  surface  node  point  for  a  given  band  point  xk  using  the  projected  distance  d ,  from 
neighbor  points  x,  to  outward  normal  line  based  on  surface  normal  vector  n  at  the 
immersed  surface  node  xs .  Large  closed  circle  represent  the  band  point  to  be  interpolated 
with  the  information  at  neighbor  point.  Hatched  black  and  gray  circles  represent  the  field 
points  and  band  points  associated  with  the  present  determination,  respectively. 


meaning  that  point  xl  is  further  away  from  the  surface  than  point  xk ,  then  the  merit  function 

12 

returns  a  very  large  value  (~H s ,  where  s  is  10"  ) 

The  actual  calculation  of  w,  is  performed  in  three  stages.  First,  only  field  points  (those 

with  d>(x;  ,t)  >  0  and  GfOCy  ,/))  =  0  )  are  considered  as  members  of  the  list  of  nearest 
neighbors.  Then,  w ,  is  calculated  according  to  Eq.  (39),  and  the  sum  of  the  weights  ^  vv/H  is 


calculated.  If  this  sum  is  non-zero,  then  the  actual  weight  function  for  each  nearest-neighbor  is 
determined  as 


(O,  = 


m 


(40) 
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Otherwise,  the  process  is  repeated,  now  considering  both  field  points  and  other  band  points  as 
members  of  the  list  of  nearest  neighbors.  If  this  application  also  results  in  no  viable 
interpolation  points  being  found,  then  the  band  point  xk  is  effectively  set  to  an  interior  point. 

The  location  at  which  interpolated  properties  are  defined,  d,  ,  is  calculated  for  a 
particular  field  point  as 

d,  =  YJ(0,(xl-xk)-h  (41) 

i 

Note  that  this  distance  is  in  the  direction  of  the  normal  coordinate.  With  this,  the  fluid  properties 
q(dj)  are  found  by  applying  the  weighting  functions 

(«) 


Flap  Struclune 


Fixed  Spar 


Upstream  Flap  Deneciiom 


Figure  4.  Aeroelastic  mesoflaps  over  a  cavity  for  (a)  subsonic  flow  and  (b)  supersonic  flow  with 

impinging  oblique  shock.  [6] 


5.  Fluid-Structure  Interactions 

One  of  the  control  devices  simulated  in  the  present  study  is  an  array  of  mesoflaps. [6-8] 
To  simulate  the  flow  control  induced  by  an  array  of  mesoflaps,  a  fluid  structure  interaction 
problem  must  solved,  as  the  mesoflaps  deform  dynamically  due  to  the  fluctuating  surface  loads. 
Figure  4  shows  a  schematic  of  a  mesoflap  system  and  its  mode  of  operation.  In  this  work,  a 
loosely  coupled  approach  is  adopted  for  the  fluid-structure  interaction  problem  with  separate 
solvers  used  for  the  fluid  and  structural  domains. 

5.1.  Structural  Solver:  Beam  Bending  Equations 

Each  mesoflap  is  treated  as  an  Euler-Bemoulli  beam  which  is  cantilevered  at  its 
upstream  end  and  deflects  due  to  the  transverse  (fluid)  pressure  loads  acting  on  its  surfaces.  The 
mesoflaps  are  rendered  as  immersed  objects  in  the  computational  domain,  and  their  motion  is 
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coupled  with  the  flow  solver  using  procedures  outlined  in  the  previous  section.  Figure  5  shows  a 
schematic  of  a  mesoflap  and  its  beam  model.  The  beam  rendered  as  1-D  mesh  (for  each  flap)  has 
101  equally  spaced  points  and  the  mesh  point  distribution  exactly  matches  that  of  the  IB  points 
on  each  mesoflap  in  the  streamwise  (x)  direction.  Two  different  governing  equations  are  used  for 
the  structural  problem,  based  on  whether  the  interaction  takes  into  account  the  inertia  of  the 
mesoflap  system  or  neglects  it.  In  general,  the  static  deflection  of  the  elastic  line  for  the  Euler- 
Bemoulli  beam  for  a  steady  load  is  given  by 


Y 


x 


Moving  part  of  mesoflap 


Fixed  Spar 


a) 


b) 


Figure  5:  a)  2-D  view  of  a  mesoflap  rendered  as  an  IB,  b)  structural  model  of  a  mesoflap  as  a 

cantilevered  beam 


d 2 

a/ 


El(x) 


d2v(x) 

d? 


f{x\  0  <x<L 


(43) 


where  v  is  the  deflection,  E  is  the  modulus  of  elasticity,  I  is  the  area  moment  of  inertia  of  the  flap 
cross-section  about  the  neutral  axis,  /  is  the  applied  transverse  force/unit  length  and  L  is  the 
beam  or  flap  length.  For  a  prismatic  beam,  the  neutral  axis  is  an  imaginary  line  passing  through 
the  centroid  of  the  cross-section  about  which  the  first  moment  of  area  (of  the  cross-section)  is 
zero.  The  beam  cross-section,  which  is  same  as  the  flap  cross-section,  is  uniform  and  rectangular 


1  3 

for  which  I  is  given  by  —  bt  ,  where  b  is  the  flap  width  and  t  is  the  flap  thickness.  As  quantities 


are  calculated  on  a  unit  width  (of  flap)  basis  in  this  work, /is  derived  on  a  per  unit  area  basis,  and 
I  is  given  by— f3 .  Since  the  beam-model  used  here  accounts  for  deformation  due  to  transverse 


loads  only,  shear  forces  on  the  mesoflap  surfaces  are  not  considered.  The  load  /  is  thus  the 
pressure  acting  on  the  flap  surface.  In  the  case  when  the  dynamics  of  the  beam  is  considered  and 
time-varying  loads  are  accounted  for,  the  governing  equation  is  that  of  a  vibrating  beam , 


dx2 


El(x ) 


d  2v{pc,t) 


dx 2 


+  f(x, t )  =  m(x )  —  —  ,0  <x<L 


dr 


(44) 


where  m[x )  is  the  mass  per  unit  length  of  the  beam.  Central  differences  are  used  for  the 
discretization  of  the  fourth  order  equation  [31]  resulting  in  a  semi-discretized  system, 

El[K’]{v}  +  pA[l]{v}  =  f  (45) 
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where  [AT']  is  the  discrete  fourth  order  differential  operator  for  the  entire  beam,  p  is  the  density 

of  the  beam  and  A  is  the  area  of  the  beam  per  unit  depth  which  is  effectively  equal  to  the  height 
of  the  beam  ( t )  in  this  case.  The  discretization  of  the  temporal  derivative  for  this  problem  is  done 
using  Newmark’s  average  acceleration  method  [32] 

v"+1  =  v"  +  At  ((1  -  y)vn  +  yvn+l ) 

1  (46) 

vn+1=vn+ A tvn  +-(1-2/?)  A  t2vn  +  j3At2vn+l 


In  the  above  equation,  the  superscript  n  denotes  the  time  step,  and  /?  and  y  are  the 
Newmark  parameters.  For  the  average  acceleration  method,  which  is  fully  implicit, 

unconditionally  stable  and  is  second  order  accurate,  /?  =  -^  ar|d  Y  =  ~  •  The  use  of  Eq.  (45)  in  Eq. 

(46)  results,  after  some  rearrangement,  in  the  following  equation, 


+ 


f-i—0 

M 

U/s  J 

l  j 

J 

(47) 


where  [jk]  =  EI[K']  is  a  coefficient  matrix  which  is  similar  to  the  stiffness  matrix  as  computed 

in  finite  element  computations.  For  the  case  of  a  quasi-steady  problem,  which  does  not  involve 
any  time  dependant  terms,  the  resulting  matrix  equation  to  solve  is  given  by, 

El[K]{vn+l}  =  f  (48) 

This  results  in  a  coefficient  matrix  which  is  not  diagonally  dominant  for  the  quasi-steady 

case  and  a  direct  solver  using  Gaussian  elimination  with  partial  pivoting  is  used  to  solve  for  the 
displacements.  For  the  dynamic  problem,  the  coefficient  matrix  is  dependent  on  the  time  step 
based  on  which  the  matrix  may  or  may  not  be  diagonally  dominant. 


5.2.  Fluid-Structure  Coupling 

5.2.1,  Determination  of  loads  on  flap  surfaces  from  neighboring  fluid  cells 

One  of  the  primary  steps  involved  in  the  coupling  of  the  fluid  and  structural  solvers  is  the 
determination  of  loads  on  the  surface  of  the  structure.  As  the  mesoflaps  are  modeled  as  Euler- 
Bemoulli  beams  in  this  study,  only  the  transverse  loads  acting  on  the  top  and  bottom  surfaces  of 
the  mesoflaps  need  to  be  determined.  The  procedure  adopted  in  determining  the  load  is  as 
follows.  The  mesoflaps  and  supporting  structures  (Figure  5)  are  rendered  as  immersed  objects  in 
the  fluid  domain  which  implies  that  there  will  be  band  points  in  the  fluid  domain  that  have  their 
nearest  surface  points  on  a  flap  surface.  So  for  each  band  point  in  the  domain,  it  is  checked  if  the 
nearest  IB  surface  point  is  a  mesoflap  upper  or  lower  surface.  If  this  is  true,  then  the  data  stored 
at  the  interpolation  point  (for  the  given  band  point)  is  projected  to  the  surface  point  on  the 
mesoflap.  In  the  case  that  a  surface  point  (on  the  mesoflap)  is  closest  to  more  than  one  band 
point,  then  the  pressure  value  from  the  interpolation  point  which  is  closer  to  the  surface  point  is 
considered.  This  procedure  results  in  a  sparsely  populated  pressure  data  set  for  each  of  the 
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mesoflap  surfaces,  which  is  then  reconstructed  using  a  bi-linear  interpolation  to  obtain  complete 
surface-load  data.  This  is  elaborated  upon  in  the  next  section. 

5.2.2.  Interpolation  of  loads 

A  bi-linear  interpolation  approach  -  linear  interpolation  with  sweeps  in  the  streamwise 
and  spanwise  direction  —  is  adopted  for  the  reconstruction  of  the  load  on  the  mesoflap  surfaces. 
This  is  done  after  a  partially  populated  load  is  obtained  on  the  mesoflap  surface  as  outlined  in  the 
previous  section.  Figure  6  shows  a  schematic  of  the  interpolation  scheme  used  to  reconstruct  the 
load  on  the  mesoflap  surfaces.  Figure  6  shows  an  array  of  IB  points  rendering  the  mesoflap 
surface,  with  ‘i’  denoting  the  streamwise  direction  and  ‘j’  the  spanwise  direction.  Load  values 
are  computed  by  first  doing  a  sweep  in  the  ‘i’  direction  and  then  in  the  ‘j’  direction.  In  the  ‘i’ 
direction,  an  interpolation  based  on  curved  surface  distances  is  used  to  compute  the  interpolated 
load  with  data  available  from  the  nearest  (populated)  right  and  left  neighbors.  This  will  populate 
the  load  values  at  all  the  points  on  those  rows  (fixed  ‘j  ’  value)  where  there  is  some  initial  data 
but  will  not  affect  the  ones  without  any  initial  load.  Thus,  following  the  ‘i’  sweep,  a  sweep  in  the 
‘j’  direction  is  performed  in  which  a  linear  interpolation  is  done  using  straight-line  distances. 
This  sweep  ensures  that  all  the  IB  points  on  the  flap  top  and  bottom  surfaces  are  populated  with 
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O  IB  point  with  load  obtained  by  interpolation 

Figure  6:  Schematic  representation  of  interpolation  scheme  used  in  the 
reconstruction  of  surface  loads  on  mesoflap  upper  and  lower  surfaces 

load  values.  The  choice  of  using  curvi-linear  distances  in  the  streamwise  (‘i’)  direction  and 
linear  distances  in  the  spanwise  (‘j’)  direction  is  made  as  the  flap  has  a  curved  shape  in  the  ‘i’ 
direction  but  is  flat  in  the  ‘j’  direction.  To  obtain  load  values  at  the  points  which  have  only  one 
populated  neighbor  (specifically  the  points  at  the  edge),  a  first  order  extrapolation  is  used. 

5.2.3.  Fluid-structure  iteration 


Once  the  interpolation  of  loads  is  done,  a  fully  populated  array  of  loads  is  obtained  for 
the  upper  and  lower  surfaces  for  each  mesoflap.  A  span-averaged  1-D  load  distribution  is  then 
obtained  for  each  mesoflap  which  is  transferred  to  the  beam  solver.  The  beam  solver  updates  the 
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deflection  and  velocity  of  the  beam  based  on  the  new  loads  by  solving  Eq.  (43)  (quasi-steady)  or 
Eq.  (44)  (dynamic),  which  is  then  used  to  obtain  a  new  configuration  for  the  mesoflap.  This 
structural  iteration  can  be  summarized  as  follows: 

•  Determination  of  loads  at  (some)  discrete  IB  points  on  the  flap  surfaces  using  band 
point  and  interpolation  point  data. 

•  Interpolation  of  loads  to  obtain  a  complete  data  set  for  the  IB  points  rendering  the 
mesoflap  upper  and  lower  surfaces. 

•  Calculation  of  a  span-averaged  1-D  load  distribution  for  the  discretized  beam. 

•  Invocation  of  the  beam  solver:  computation  of  new  position  and  velocities  of  mesoflap 

•  Reclassification  of  fluid  domain  into  field,  band  and  interior  cells. 

The  frequency  at  which  a  structural  iteration,  is  done  depends  on  whether  a  quasi-steady 
or  fully  dynamic  fluid-structure  interaction  problem  is  being  solved.  For  the  quasi-steady  2-D 
simulation  this  is  done  after  every  2500  iterations  of  the  fluid  solver  which  is  considered 
sufficient  to  obtain  a  reasonably  converged  intermediate  fluid  solution.  The  process  of  updating 
the  configuration  of  the  mesoflaps  (after  every  2500  fluid-solver  cycles)  is  continued  until  a 
convergence  criterion  is  met  for  the  computed  displacements.  This  is  evaluated  in  the  following 
way: 

•  At  every  structural  iteration,  the  maximum  change  in  displacement  of  each  mesoflap 
(compared  to  the  previous  configuration)  is  computed. 

•  It  is  then  checked  whether  this  value  is  less  than  a  specified  tolerance,  which  is  chosen 
as  10%  of  the  flap  thickness  [32],  If  this  is  true  for  all  the  mesoflaps  then  the  fluid- 
structure  interaction  is  assumed  to  be  converged 

•  A  final  fluid  solution  is  obtained  for  the  converged  configuration  of  the  mesoflap 
system. 

For  the  fully  dynamic  approach,  a  structural  iteration  is  done  at  each  time  step  (every  fluid 
iteration). 

5.2.4.  Fluid-solver  changes 

A  couple  of  changes  are  adopted  for  the  fluid  solver  to  facilitate  the  smooth  running  of 
the  fluid-structure  interaction  problem  when  the  dynamics  of  the  structure  is  considered.  These 
are  listed  below 

•  A  reduced  order  of  time  integration  (Euler  implicit)  is  used  for  the  fluid  cells  in  the 
vicinity  of  the  immersed  surfaces. 

•  For  the  baseline  IB  method  as  applied  to  non-moving  geometries,  the  solution  in  the 
interior  cells  is  fixed  somewhat  arbitrarily  and  not  as  function  of  the  solution  of 
neighboring  cells.  This  approach  was  making  the  solution  unstable  for  the  moving  IB 
surfaces  investigated  herein.  A  possible  reason  for  this  is  as  follows.  As  the  continuity 
equation  is  integrated  for  the  band  cells  in  this  method,  this  requires  the  knowledge  of 
flow  properties  stored  at  these  band  cells  at  the  previous  time  step  (for  updating  and 
residual  construction).  For  the  case  in  which  a  cell  changes  from  being  an  interior  cell  to  a 
band  cell,  the  values  being  used  for  the  residual  construction  may  not  be  consistent.  To 
address  this  issue,  the  reconstruction  of  the  solution  in  the  interior  cells  which  are 
neighbored  by  a  band  cell  has  been  modified.  If  any  interior  cell  has  as  its  neighbor  a 
band  cell,  then  its  properties  are  evaluated  by  interpolating  the  values  of  its  neighboring 
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band  cells  using  an  inverse  distance  weighting  approach.  This  is  similar  in  concept  to  the 
treatment  used  for  simulation  of  moving  bodies  by  Balaras  [33]. 


6.  Results 

6.1.  Overview 

Three  main  sets  of  experiments  are  simulated  in  this  work.  The  first  set  of  experiments 
was  conducted  at  Cambridge  University  [10]  and  involved  supersonic  flow  over  single  micro 
VGs  as  well  as  impinging  oblique  shock  interactions  with  and  without  micro  VG  effects.  The 
second  set  of  experiments  was  conducted  at  NASA  Glenn  Research  Center  [26,34]  and  involved 
flat  plate  boundary  layer  flows  with  bleed  effects  and  impinging  oblique  shock  interactions  with 
and  without  bleed.  The  third  set  of  experiments  was  performed  at  the  University  of  Illinois  [6-8] 
and  involved  impinging  oblique-shock  interactions  with  and  without  meso-flap  flow  control. 
Results  from  simulations  of  each  of  these  experiments  are  described  in  detail  in  the  following 
sections.  A  collaboration  with  NASA  Langley  led  to  the  use  of  the  immersed-boundary  flow 
solver  to  investigate  linear  stability  of  laminar  flows  influenced  by  boundary  layer  trip  elements. 
This  work  will  be  discussed  in  brief,  as  will  results  from  a  workshop  study  that  applied  the  IB 
method  to  simulate  the  effects  of  a  shock  generator. 

Free-stream  conditions  and  inflow  boundary  layer  properties  for  the  three  main  sets  of 
experiments  are  summarized  in  Table  2. 


Tal 

3le2:  Inflow  boundary  layer  pro] 

jerties 

Parameter 

Babinsky,  et  al. 

Willis,  et  al. 

Gefroh,  et  al. 

2.5 

2.46 

2.41 

Po  (Pa) 

3.056e6 

1.724e6 

5.026e6 

T0  (K) 

290.0 

293 

300 

Re/m 

30.0e6 

17.5e6 

50.4e6 

S0  (cm) 

0.6 

2.63 

0.4 

6.2.  Micro  VG  Control  of  Shock  /  Boundary  Layer  Interactions 

The  first  set  of  calculations  described  in  this  work  correspond  to  experiments  that 
examine  the  effects  of  micro  vortex  generators  (micro  VGs)  in  controlling  shock  /  boundary 
layer  interactions.  Here,  the  individual  vortex  generators  (up  to  3)  are  rendered  as  immersed 
objects.  The  following  sections  are  excerpted  from  papers  P-1  and  P-5  in  the  publication  list 
(Section  9). 
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6.2.1.  Experiment  details 


Experiments  involving  the  effects  of  micro  VGs  on  shock  interaction  control  have  been 
performed  at  Cambridge  University  by  Dr.  Holger  Babinsky  and  his  students.  The  experiments 
were  performed  at  a  nominal  Mach  number  of  2.5,  a  stagnation  temperature  of  290  K,  and  a 
Reynolds  number  /  meter  of30  x  106 .  The  test  section  in  the  wind  tunnel  is  90  mm  high  and  1 10 
mm  wide.  Data  collected  in  these  experiments  includes  pitot  pressure  surveys,  wall  static 


Figure  7:  Schematic  of  Cambridge  University  wind  tunnel  showing  the  locations  of  the  data 
stations,  VG  positions,  and  shock  generator  (provided  by  H.  Babinsky) 

pressure  distributions,  and  axial  velocity  profiles  obtained  at  various  streamwise  and  spanwise 
stations  using  laser  Doppler  anemometry  (LDA).  Of  the  cases  considered  in  the  Cambridge 
database,  this  paper  focuses  on  three  experiments: 

1.  Mach  2.5  flow  over  single  wedge-shaped  micro  VGs  of  heights  h  =  3  mm  and  h  =  6  mm 
(LDA  velocity  data,  Schlieren  imaging,  surface  oil-flow) 

2.  An  oblique  shock  /  turbulent  boundary  layer  interaction  at  Mach  2.5  induced  by  a  7- 
degree  wedge  placed  on  the  upper  wall  of  the  wind  tunnel  (LDA  velocity  data,  surface 
pressure  measurements,  Schlieren  imaging,  surface  oil-flow) 

3.  The  same  oblique  shock  /  turbulent  boundary  layer  interaction  with  the  addition  of  an 
array  of  three  micro  VGs  with  heights  h  =  3  mm  upstream  of  the  shock  impingement  point 
(LDA  velocity  data,  surface  pressure  measurements,  Schlieren  imaging,  surface  oil-flow) 


24 


Figure  7  shows  the  general  arrangement  of  the  VG  arrays  within  the  Cambridge  wind 
tunnel.  For  the  experiments  considered  herein,  the  trailing  edge  of  each  micro  VG  was  located 
108.5  mm  downstream  of  the  start  of  the  test  section.  The  leading  edge  of  the  shock  generator 
was  also  placed  at  the  start  of  the  test  section,  resulting  in  a  lower-wall  shock-impingement 
position  of  159  mm,  assuming  inviscid  flow.  Velocity  profdes  were  collected  at  streamwise 
stations  of  X  =  20  mm,  X  =  80  mm,  and  X  =  140  mm,  with  the  X  coordinate  measured  from  the 
trailing  edge  of  the  micro  VG.  Additional  profiles  were  collected  at  X  =  -20  mm  (in  the  absence 
of  the  VG)  to  determine  the  inflow  boundary  layer  characteristics,  and  at  X  =  50  mm  for  the 
shock-impingement  cases.  The  nominal  boundary  layer  thickness  at  the  inflow  plane  (X  =  -20 
mm)  is  around  6.5  mm.  At  each  of  the  streamwise  stations,  LDA  data  was  recorded  at  three  or 
four  spanwise  stations  (with  Z  =  0  mm  being  the  tunnel  centerline)  as  indicated  in  Table  3.  The 
dimensions  and  alignment  of  the  wedge-shaped  micro  VGs  in  the  experiments  are  based  on  [9] 
and  [10]. 


Table  3.  Spanwise  location  of  data  stations  for  different  micro  VG  sizes  with  /  without  shock 

interaction 


Experiment 

Spanwise  data  locations  (Z,  mm) 

Mach  2.5  flow  over  6  mm  micro  VG 

0.0,  4.5,  9.0,  18.0 

Mach  2.5  flow  over  3mm  micro  VG 

0.0,  2.3,  4.5,  9.0 

Mach  2.5  shock  /  boundary  layer  interaction  (3x3  mm 
micro  VG  array) 

0.0,  4.5,  9.0 

6.2.2.  Computational  domain  and  calculation  details 

Two  computational  domains  have  been  used  for  the  simulations  described  in  this  work. 
The  first  domain,  used  for  the  single  micro-VG  calculations  and  for  the  ‘idealized’  SBLI 
simulations  described  later,  ranges  from  X  =  -108.5  mm  to  X  =262.75  mm  in  the  streamwise 
direction  (with  the  X  coordinate  measured  with  respect  to  the  start  of  the  wind-tunnel  test 
section),  from  Y  =  0  mm  to  Y  =  90  mm  in  the  normal  direction,  and  from  Z  =  -22.5  mm  to  Z  = 
22.5  mm  in  the  spanwise  direction.  The  mesh  spacing  in  the  X  and  Z  directions  is  0.5  mm,  and 
a  minimum  spacing  of  0.005  mm  is  enforced  at  the  lower  wall.  The  total  mesh  size  is 
747  x  200  x  90  =  13.45  million  interior  mesh  cells.  Slip  wall  conditions  are  prescribed  on  the 
upper  boundary  (Y  =  90  mm),  extrapolation  conditions  are  prescribed  at  the  outer  boundary  ( X  = 
262.75  mm),  periodic  boundary  conditions  are  applied  at  Z  =  +/-  22.5  mm,  and  no-slip, 
adiabatic- wall  boundary  conditions  are  applied  at  Y  =  0  mm.  For  the  SBLI  simulations,  the 
shock-generator  leading  edge  is  placed  at  V  =-13  mm,  leading  to  a  lower- wall  impingement 
position  of  X  =  146  mm.  This  value  is  consistent  with  Schlieren  imaging  data  (H.  Babinsky, 
personal  communication)  and  with  full-wind  tunnel  simulations  (discussed  later)  and  accounts 
for  the  effects  of  the  upper-surface  boundary  layer  in  modifying  the  shock  inception  point. 

The  second  computational  domain  renders  the  Cambridge  blow-down  wind  tunnel  in  two 
parts.  The  first  part  encompasses  one  half  of  the  entire  wind  tunnel,  including  the  nozzle  and 
test  section.  Coordinates  of  the  nozzle  were  provided  by  Dr.  Babinsky  and  were  spline-fitted  for 
use  in  the  mesh-generation  exercise.  The  mesh  is  clustered  to  the  upper  wall,  the  lower  wall  (Y  = 
0  mm),  and  to  one  side  wall  (Z  =  55  mm)  (centerline  symmetry  being  assumed)  with  a  minimum 
spacing  of  0.005  mm.  Uniform  mesh  spacing  at  a  value  of  0.5  mm  is  used  from  Z  =  0  to  Z  = 
~50  mm,  and  a  power-law  mesh-clustering  function  is  then  used  to  connect  this  spacing  to  the 
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wall  spacing  of  0.005  mm.  The  mesh  spacing  in  the  X  direction  is  1.0  mm,  and  the  total  number 
of  interior  mesh  cells  is  840  x  225  x  130  =  24.57  million  cells.  A  steady  RANS  calculation  was 
performed  on  this  mesh  to  obtain  accurate  inflow  conditions  for  refined  simulations  of  the  flow 
in  the  test  section.  The  test-section  computational  domain  extends  from  X  =  -121.75  mm  to  J  = 
252.75  mm  and  has  a  spacing  of  0.5  mm  in  the  streamwise  direction.  The  7-degree  shock 
generator  is  placed  on  the  upper  wall  in  this  configuration  (with  the  leading  edge  at  the  start  of 
the  test  section  as  in  the  experiment).  The  numbers  of  mesh  points  and  the  stretching  patterns  in 
the  Y  and  Z  directions  are  identical  to  that  used  in  the  full  wind-tunnel  grid,  and  the  total  number 
of  interior  mesh  cells  is  750  x  225  x  130  =  21.94  million  cells.  The  wedge-shaped  micro  VGs  are 
rendered  separately  as  immersed  boundaries  using  a  simple  structured  surface  mesh  for  each  VG. 
The  present  implementation  of  the  distance-function  calculation  requires  a  relatively  fine  surface 
mesh  for  accurate  results.  For  the  VG  array  (three  3  mm  VGs),  a  total  of  2234907  surface  points 
is  used,  while  for  the  6  mm  VG  array,  a  total  of  2964240  points  is  used. 

Convergence  of  the  RANS  calculations  was  ascertained  using  three  measures:  relative 
decrease  in  residual  norm,  constancy  of  surface  quantities,  and  constancy  of  the  global  mass  flow 
rate.  All  RANS  cases  reached  a  steady-state  solution.  For  the  cases  performed  on  the  idealized 
domain,  a  flat-plate  simulation  at  the  test-section  conditions  was  performed  to  determine  the 
place  at  which  the  predicted  boundary  layer  properties  most  closely  matched  the  experimental 
values  at  X  =  -20  mm.  The  inflow  plane  used  in  the  calculations  was  then  extracted  from  the 
flat  plate  solution  at  a  location  88.5  mm  upstream  of  the  matching  position.  The  LES/RANS 
calculations  were  initiated  by  super-imposing  scaled  boundary-layer  fluctuations  from  an  earlier 
calculation  onto  part  of  a  converged  RANS  solution.  After  several  flow-through  times  to 
eliminate  initial  transients,  time-averaged  statistics  for  the  LES/RANS  calculations  were 
collected  over  a  minimum  of  0.007s  (10.7  flow-through  times  based  on  a  domain  length  of 
371.25  mm  and  a  ffee-stream  velocity  of  570  m/s). 

The  different  cases  considered  in  this  part  of  the  study  are  summarized  in  Table  4. 


Table  4:  Summary  of  cases:  SBLI  with  micro  VGs 


Case 

Domain 

Grid  size 

Turbulence 

Model(s) 

IB  details 

1.  Mach  2.5  flow  over  6  mm 
micro  VG 

Periodic 

13.45  M 

RANS  (Menter 
BSL), 

LES/RANS 

£=1,1/7, 1/9 
(compressible) 

£=1/7  (compressible, 
mass  conservative) 

2.  Mach  2.5  flow  over  3  mm 
micro  VG 

Periodic 

13.45  M 

RANS  (Menter 
BSL), 

LES/RANS 

£=1/7  (compressible, 

mass 

conservative) 

3.  Mach  2.5  SBLI 

Wind 

Tunnel 

21.94  M 

RANS  (Menter 
SST) 

None 

4.  Mach  2.5  SBLI  with  micro 
VG  array 

Wind 

Tunnel 

21.94  M 

RANS  (Menter 
SST) 

£=1/7  (compressible, 

mass 

conservative) 

5.  Mach  2.5  SBLI 

Periodic 

13.45  M 

RANS  (Menter 
SST), 

LES/RANS 

None 

6.  Mach  2.5  SBLI  with  micro 
VG  array 

Periodic 

13.45  M 

RANS  (Menter 
SST), 

LES/RANS 

£=1/7  (compressible, 

mass 

conservative) 
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6.2.3.  Inflow  boundary  layer 


Inflow  velocity  profiles  were  measured  at  X  =  -20  mm  in  the  absence  of  the  micro  VG 
array.  Figure  8  compares  inflow  velocity  profiles  obtained  on  the  idealized  domain  (RANS  and 
LES/RANS),  the  wind  tunnel  domain 


(RANS),  and  in  the  experiment. 
Compared  with  the  experiment,  the 
computed  profiles  show  consistently 
larger  values  of  velocity  in  the  inner  part 
of  the  boundary  layer,  which  implies 
lower  values  of  the  displacement  and 
momentum  thicknesses.  The  velocity 
profile  obtained  from  the  wind-tunnel 
calculation  is  somewhat  thicker  than  the 
profiles  obtained  on  the  idealized  domain. 
The  LES/RANS  and  RANS  predictions 
are  very  similar,  with  the  LES/RANS 
profile  showing  slightly  higher  values  of 
velocity  in  the  near-wall  region. 


6.2.4.  Immersed-boundary  model 
assessment:  Mach  2.5  flow  over  single 
micro  VG 


U  (m/s) 

Figure  8:  Inflow  boundary  layer  profiles 
(X  =  -20  mm) 


To  assess  the  immersed  boundary  technique,  results  obtained  using  variations  of  the  IB 
method  are  compared  with  a  solution  obtained  on  a  body  fitted  grid  and  with  experimental  data 
for  Mach  2.5  flow  over  a  6  mm  micro  VG.  The  body-fitted  grid  was  generated  using  the 
commercial  software  GRIDGEN  and  contains  13.2  million  cells,  clustered  to  each  solid  surface 
with  a  minimum  spacing  of  0.005  mm.  The  X  extent  of  the  body- fitted  mesh  ranges  from  -108.5 
mm  to  224.5  mm.  Contour  plots  of  Mach  number  at  X  =  2  mm  and  X  =  90  mm  and  centerline 
plots  of  axial  velocity  at  X  =  20  mm  and  X  =  80  mm  are  used  in  the  assessment.  The  scale  for 
the  contour  plots  ranges  from  0  to  2.85,  with  0  being  the  darkest  contours. 

Mach  number  contours  just  downstream  of  the  micro  VG  trailing  edge  at  X  =  2  mm 
(Figure  9)  clearly  show  the  initial  growth  of  the  primary  vortex  cores.  The  spanwise  extent  of  the 
low  momentum  region  is  more  pronounced  for  the  k=  1  interpolation  and  is  least  for  the  solution 
obtained  on  the  body-fitted  mesh.  The  vortices  lift  off  less  in  the  IB  solutions  and  spread  more 
laterally  compared  to  the  body-fitted  mesh  solution.  Also,  the  higher-momentum  flow  reaches 
closer  to  the  wall  for  the  body-fitted  mesh  solution.  This  implies  a  stronger  interaction  between 
the  vortex  cores  and  the  external  flow.  The  presence  of  a  secondary  pair  of  vortices  formed  in 
the  near- wall  region  downstream  of  the  micro  VG  trailing  edge  is  indicated  by  darker  contours 
near  the  surface.  The  exact  shape  of  the  wake  as  predicted  on  the  body-fitted  grid  is  not  very  well 
matched  by  any  of  the  RANS  IB  solutions,  but  a  trend  toward  improved  results  with  lower 
values  of  k  is  evident.  LES/RANS  results  using  k=\/l  and  the  mass-conservative  IB  are  the 
closest  to  the  body-fitted  mesh  RANS  solutions. 

At  X  =  90  mm  (Figure  10),  the  vortex  pairs  generated  behind  the  micro  VG  spread  out 
and  lift  from  the  surface,  and  the  lower-momentum  core  region  becomes  less  pronounced,  as 
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Figure  9:  Mach  contours  at  X  =  2  mm  for  6  mm  VG  placed  in  supersonic 

boundary  layer 


suggested  by  the  lighter  contours.  The  body-fitted  solution  again  shows  a  narrower  and  taller 
low-momentum  region,  but  the  differences  between  this  solution  and  the  IB  solutions  are  less 
apparent  than  at  X  =  2  mm.  The  low-momentum  region  is  larger  for  the  Ic=  1  interpolation  than 
for  the  others,  implying  that  this  ‘laminar’  choice  might  not  provide  enough  flow  attachment  to 
the  micro  VG  surface.  The  LES/RANS  solution  shows  generally  higher  velocities  in  the  wake 
and  close  to  the  lower  wall,  which  is  indicative  of  a  stronger  vortex  providing  more  entrainment 
of  high-momentum  fluid. 

Centerline  velocity  profiles  in  Figure  1 1  show  that  the  body- fitted  grid  solution  seems  to 
match  the  experimental  data  best  for  both  the  stations  shown.  The  IB  predictions  approach  the 
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body-fitted  solution  for  lower  values  of  k,  and  interestingly,  the  mass-conservative  IB  method  is 
not  as  quite  as  accurate  as  the  others.  From  the  comparisons  at  the  Z  =  20  mm  station,  it  appears 
that  the  effects  of  the  primary  vortex  in  energizing  the  near-wall  boundary  layer  are  somewhat 
under-predicted  by  all  of  the  IB  methods.  The  shape  of  the  wake  deficit  is  not  predicted  well  by 
any  of  the  RANS  methods  at  X  =  20  mm  but  improves  somewhat  at  X  =  80  mm.  These  results 


z  (mm) 


z  (mm) 


Figure  10:  Mach  contours  at  x  =  90  mm  for  6  mm  VG  in 
supersonic  boundary  layer 


indicate  that  the  IB  rendering  of  the  micro  VG  provides  more  of  a  momentum  ‘sink’  than  does 
the  body-fitted  rendering.  Considering  that  no  effort  is  made  to  resolve  the  flow  near  the 
surfaces  of  the  micro  VG  when  the  IB  techniques  are  used,  the  level  of  agreement  shown  is  still 
encouraging  enough  for  further  evaluation. 
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Figure  11:  Comparison  of  centerline  axial  velocity  profiles  for  different 
formulations  of  the  IB  method  (Menter’s  k-co),  a  body  fitted  grid  (Menter’s  k-  co) 

and  experimental  data. 

6.2.5.  Detailed  data  comparisons  for  Mach  2.5  flow  over  6mm  and  3mm  micro  VGs 

Figure  12  compares  velocity  profiles  obtained  using  the  RANS  and  LES/RANS  models 
(k=\H,  mass-conservative  IB)  with  LDA  data  collected  at  several  X  and  Z  positions  behind  the 
6mm  micro  VG.  (see  Table  1)  The  energizing  effects  of  the  primary  vortices  are  better  captured 
by  the  LES/RANS  model  than  by  the  RANS  model.  This  is  reflected  in  the  fuller  velocity 
profiles  near  the  wall  and  close  to  the  VG  (Z  =  0  mm  and  Z=  4.5  mm),  which  agree  quite  well 
with  the  experimental  data.  However  there  is  a  considerable  disagreement  in  the  level  of  the 
wake  deficit  found  at  the  Z  =  4.5  mm  station  at  X  =  20  mm  and  X  =  80  mm.  In  the  experiment, 
the  wake  deficit  is  almost  absent  at  this  Z  station,  either  implying  less  lateral  spreading  of  the 
vortices  or  a  device  misalignment  that  shifts  the  time-mean  position  of  the  vortex  pair.  It  should 
be  noted  that  all  measurements  were  taken  on  one  side  of  the  centerline.  In  an  attempt  to 
determine  whether  uncertainties  in  the  placement  of  the  laser  might  play  a  role,  the  LES/RANS 
velocity  profiles  were  also  averaged  in  the  Z  direction  over  the  experimental  positional 
uncertainty  of  2.5  mm  (H.  Babinsky,  personal  communication).  These  profiles  show  a  slight 
improvement  with  respect  to  the  experimental  data.  The  good  agreement  shown  at  the  most 
downstream  station  (X  =  140  mm)  for  the  LES/RANS  model  lends  support  to  the  possibility  of  a 
slight  misalignment  of  the  micro  VG  with  respect  to  the  incoming  flow.  The  level  of  agreement 
found  at  the  more  outboard  stations  (Z  =  9  mm  and  Z  =  1 8  mm)  is  generally  good,  though  near¬ 
wall  discrepancies  similar  to  those  observed  for  the  inflow  velocity  profile  persist.  The 
LES/RANS  predictions  are  in  better  agreement  with  experimental  data  at  most  stations. 

Figure  13  presents  a  similar  set  of  comparisons  for  Mach  2.5  flow  over  a  3  mm  micro 
VG.  Centerline  velocity  profiles  from  the  RANS  and  LES/RANS  models  do  not  agree 
particularly  well  with  the  data,  in  that  the  latter  shows  a  much  smaller  wake  deficit.  The 
spanwise  extent  of  the  vortex  pair  is  predicted  to  be  larger  than  indicated  in  the  experiment,  as 
the  measured  wake  deficit  seems  not  to  exist  at  any  station  other  than  the  centerline.  Agreement 
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with  experimental  data  further  away  from  the  vortex  pair  is  satisfactory  for  both  models  and 
again,  the  LES/RANS  model  provides  more  accurate  predictions  at  most  locations. 
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Figure  12:  Comparison  of  axial 
velocity  profiles  for  flow  over  6 
mm  VG  in  idealized  domain  (*  = 
data  averaged  over  a  span  wise 
filter  of  2.5  mm,  centered  at  the  Z 
location) 


Figure  13:  Comparison  of  axial 
velocity  profiles  for  flow  over  3 
mm  VG  in  idealized  domain 
(*  =  data  averaged  over  a  span 
wise  filter  of  2.5  mm,  centered  at 
the  Z  location) 
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6.2.6.  Wind-tunnel  shock  /  boundary  layer  interaction  -  flow  structure 


Figures  14  and  15  present  near-surface  axial  velocity  contours  and  streamtraces  for  the 
wind-tunnel  shock  /  boundary  layer  interaction  without  micro- VG  flow  control.  The  vertical 
white  lines  in  the  velocity  contour  plots  indicate  the  theoretical  position  of  the  shock 
impingement  (X  =  159  mm),  while  the  zero-velocity  contour  is  indicated  by  a  thinner  white  line. 
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Figure  14:  Near  surface  velocity  contours  at  Y  =  0.0025  mm  (range:  -5 
m/s  (dark)  to  35  m/s)  for  SBLI  in  wind  tunnel  (contours  reflected  about 

centerline  for  clarity) 

The  interaction  of  the  oblique  shock  with  the  sidewall  boundary  layer  results  in  a  crossflow 
separation  similar  to  that  observed  in  sharp  fin  interactions.  The  downward  motion  of  the 
crossflow  is  arrested  by  the  flat  plate,  leading  to  the  formation  of  vortical  structures  located  in 
the  sidewall  /  flat  plate  junctures  (indicated  in  dark  contours  near  the  upper  and  lower  walls). 
The  displacement  effects  of  these  structures  induce  the  formation  of  a  pair  of  separation  shocks 
and  a  general  movement  of  the  flat-plate  boundary  layer  fluid  toward  the  centerline.  Part  of  the 
near-wall  fluid  lifts  away  from  the  surface  at  two  foci  of  separation,  located  at  the  points  of 
intersection  of  the  separation  shocks  with  the  impinging  oblique  shock.  (Figure  15)  This  fluid 
then  spirals  through  the  main  separation  region  before  exiting  near  the  centerline.  The  growth 
of  the  separation  region  toward  the  sidewall  is  prevented  by  the  crossflow  induced  by  the  comer 
vortices,  leading  to  a  non-uniform  region  of  low-momentum  fluid  that  expands  more  near  the 
centerline.  (Figure  14) 

6.2.7.  Wind-tunnel  shock  /  boundary  layer  interaction  -comparisons  with  experimental  data 
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Figure  15  compares  centerline  surface  pressure  distributions  for  the  SBLIs  with  and 
without  micro  VG  control  with  experimental  data.  The  separation  shock  strength  and  the  level 
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of  upstream  influence  of  the  interaction  are  predicted  well  by  the  Menter  SST  model  for  the  case 
without  the  micro  VG  array.  The  predicted  pressure  distribution  shows  the  start  of  a  plateau 
behind  the  initial  pressure  rise.  This  feature  is  commonly  observed  in  stronger  SBLIs  but  is  not 
seen  in  the  experimental  distribution.  Velocity  profiles  at  different  streamwise  stations  along  the 
centerline  are  shown  in  Figure  17. 

The  first  station  (X  =  20  mm)  is 
within  the  region  of  upstream 
influence,  and  the  predicted  velocity 
profile  shows  a  small  region  of 
reversed  flow  that  is  not  captured  by 
the  LDA  data.  The  comparison  at 
X  =  50  mm  indicates  that  the 
calculation  does  not  quite  capture 
the  shape  of  the  wake-like  velocity 
profile  near  the  re-attachment 
position.  The  recovery  of  the 
centerline  velocity  (X  =  80  mm  and 
X  =  140  mm)  is  predicted  very  well, 
except  for  some  deviations  in  the 
free-stream  velocity  similar  to  those 
noted  earlier. 
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6.2.8.  Wind-tunnel  shock  /  boundary 

layer  interaction  with  micro- VG 

flow  control  -  flow  structure 


Figure  16:  Centerline  surface  pressure  distributions  for 
wind-tunnel  cases 


The  effects  of  the  micro  VG  array  on  the  near-surface  flow  structure  are  illustrated  in 
Figure  18  (surface  velocity  contours)  and  in  Figure  19  (near- wall  streamtraces).  Darker  regions 
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in  Figure  1 8  represent  slower-moving  fluid,  while  lighter  regions  represent  faster-moving  fluid. 


MenterSST 

Experiment 


U(nVs) 

Figure  17:  Centerline  axial  velocity  profiles 


wind-tunnel  case  without  micro  VG  array 


The  major  effect  of  the  micro  VG  array  is  to  induce  the  formation  of  longitudinal  vortical 
structures  (one  pair  per  micro  VG).  These  structures  force  higher  momentum  fluid  toward  the 
surface,  energizing  the  inner  part  of  the  boundary  layer.  The  traces  of  the  vortical  structures  are 
indicated  as  light  bands  in  Figure  18.  The  region  directly  behind  the  trailing  edge  of  each  micro 
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Figure  18:  Near-surface  velocity  contours  at  Y  =  0.0025  mm  (range:  -5  m/s  (dark)  to  35 
m/s)  for  SBLI  with  micro  VG  array  in  wind  tunnel  (contours  reflected  about  the 

centerline  for  clarity) 

VG  is  influenced  by  another  pair  of  counter-rotating  vortices  that  act  to  move  fluid  away  from 
the  surface,  resulting  in  low-momentum  streaks  that  also  persist  downstream  of  the  interaction 


streaks  of  high-momentum  fluid 


streaks  of  lew-momentum  fluid  behind  VG 
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region.  Comparing  Figures  18  and  14,  it  is  clear  that  the  vortices  induced  by  the  micro  VG 
array  deform  the  large  separation  region  but  do  not  significantly  reduce  its  size.  The  level  of 
upstream  influence  is  reduced  in  regions  where  the  counter-rotating  vortices  energize  the  near¬ 
wall  boundary  layer  but  increases  in  the  regions  between  the  micro  VGs.  The  strengths  of  the 
vortices  generated  by  the  centerline  micro  VG  appear  to  be  reduced  relative  to  the  others,  as  a 
consequence  of  their  interaction  with  the  large  separation  region.  The  amount  of  reversed  flow 
again  approximately  coincides  with  the  outer  boundary  of  the  dark  region,  except  for  the  portions 
directly  influenced  by  the  counter-rotating  vortex  pairs. 

The  black  streamtraces  in  Figure  19  show  that  part  of  the  near- wall  fluid  originating 
between  the  outer  VGs  and  the  (corresponding)  sidewalls  is  entrained  into  the  vortex  pairs 
generated  from  these  devices.  Another  portion  moves  toward  the  centerline  and  lifts  from  the 
surface  at  two  foci  of  separation,  where  it  is  entrained  into  the  vortices  generated  by  the 
centerline  micro  VG.  Some  of  the  near- wall  fluid  originating  between  the  micro  VGs  is  also 


Figure  19:  Near-surface  streamtraces  for  SBLI  with  micro  VG  array 

in  wind  tunnel 

entrained  into  the  low-momentum  region  and  departs  from  the  surface  at  these  foci.  Fluid 
originating  further  away  from  the  surface  (red  streamtraces)  moves  over  and  around  the  micro 
VGs  but  does  not  enter  the  low-momentum  region.  Some  of  this  fluid  eventually  migrates  into 
the  vortical  structures. 

6.2.9.  Wind-tunnel  shock  /  boundary  layer  interaction  with  micro-VG  flow  control  - 

comparisons  with  experimental  data 

Pressure  distributions  shown  in  Figure  16  indicate  that  the  effect  of  the  micro  VG  array, 
at  least  at  the  centerline,  is  to  reduce  slightly  the  upstream  extent  of  the  separation  region  and  to 
sharpen  the  pressure  rise.  An  overshoot  in  pressure  is  observed  near  the  re-attachment  location 
for  both  the  calculated  and  measured  distributions  but  is  more  pronounced  in  the  former. 

Axial  velocity  predictions  at  several  streamwise  and  spanwise  stations  are  shown  in 
Figure  20.  At  the  X  =  20  mm  and  X  =  50  mm  stations  along  the  centerline  (Z  =  0),  the 
calculated  velocity  profiles  show  a  pronounced  low-momentum  region  that  is  associated  with  the 
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Figure  20:  Axial  velocity  profiles  for  wind-tunnel  case  with  micro  VG  array  (*  indicates  data 
averaged  over  a  spanwise  filter  of  2.5  mm,  centered  at  the  Z  location) 


vortex  pairs.  This  feature  is  not  prominent  in  the  experimental  profiles,  similar  to  the  single 
micro-VG  results.  The  centerline  profiles  further  downstream  are  in  better  agreement  with  the 
experimental  data.  The  Z  =  4.5  mm,  X  =  50  mm  position  is  near  the  location  where  the 
centerline  vortex  divides  the  main  separation  region  into  a  smaller  one  located  near  the  centerline 
and  a  larger  one  located  further  away.  The  topology  of  the  flow  in  this  region  changes  abruptly, 
and  improved  predictions  are  obtained  by  averaging  the  computed  velocity  profiles  over  the 
experimental  positional  uncertainty  of  2.5  mm.  The  same  is  true  for  the  data  station  at  Z  =  9.0 
mm,  X  =  50  mm.  Agreement  with  the  experimental  data  at  X  =  80  mm  and  X  =  140  mm  can  be 
considered  as  excellent. 


6.2.10.  Idealized  shock  /  boundary  layer  interaction 


It  is  clear  from  the  preceding  discussion  that  three-dimensional  effects  significantly 
influence  the  flow  features  in  the  Cambridge  wind  tunnel  experiments,  and  the  evaluation  of  the 
expected  benefits  of  the  micro  VG  array  is  made  more  complicated  as  a  result.  To  assess  the 
behavior  of  the  micro  VG  concept  in  an  idealized  setting  and  to  compare  the  predictions  obtained 
from  the  RANS  and  LES/RANS  models,  we  consider  an  oblique  shock  interaction  with  and 
without  micro  VG  control  using  the  periodic  domain  discussed  earlier.  Here,  the  viscous  side 
and  top  walls  of  the  wind  tunnel  are  neglected,  and  the  incoming  boundary  layer  is  developed  as 
a  flat-plate  solution  at  the  test-section  conditions.  Periodic  boundary  conditions  are  imposed  in 
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the  spanwise  (Z)  direction.  The  shock  generator  position  is  adjusted  so  that  the  oblique  shock 
impinges  at  X  =  -146  mm  as  discussed  earlier.  The  resolution  of  the  interaction  region  in  the 
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Figure  21:  Centerline  surface  pressure 
distributions  for  idealized  cases 

streamwise  and  spanwise  directions  is  the  same  as  in  the  wind-tunnel  calculations.  Both  RANS 
and  hybrid  LES/RANS  calculations  are  performed  for  this  configuration. 

In  the  absence  of  the  sidewalls  and  the  VG  array,  the  oblique  shock  interaction  is 
nominally  two-dimensional.  Centerline  pressure  distributions  shown  in  Figure  21  imply  that  the 
upstream  extent  of  the  low-momentum  region  is  much  less  in  the  idealized  interaction  than  in  the 
experiment.  As  such,  one  cannot  expect  quantitative  agreement  with  the  experimental  velocity 
profiles  in  the  interaction  region,  and  no  such  comparisons  are  presented.  Details  of  the 
idealized  interaction  are  shown  in  the  near-surface  velocity  contours  of  Figures  22  and  23.  To 
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Figure  22:  Surface  velocity  contours  (range:  -5  m/s  (dark)  to  35  m/s)  for 
idealized  SBLI  (Menter  SST  RANS  model) 
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facilitate  a  direct  visual  comparison  with  the  wind-tunnel  interaction,  the  Z  scale  is  expanded  to 
approximate  the  actual  width  of  the  wind  tunnel,  and  both  the  velocity  contour  levels  and  the  X 
scale  are  the  same  as  in  Figure  14.  The  extent  of  the  low-momentum  region  is  captured  similarly 
by  both  modeling  approaches,  but  the  LES/RANS  model  predicts  a  smaller  region  of  reversed 
flow  with  smaller  velocity  components  (in  magnitude).  The  near-surface  boundary  layer 
recovers  more  rapidly  in  the  LES/RANS  calculation,  reaching  span-averaged  values  of  ~27  m/s 
as  compared  with  ~20  m/s  for  the  RANS  solution  at  the  exit  plane. 

One  obvious  feature  is  that  traces  of  organized  three-dimensionality  are  indicated  in  the 
time-averaged  LES/RANS  contours. (Figure  23)  These  result  from  the  combination  of  several 
factors,  the  most  noteworthy  of  which  is  the  tendency  of  the  recycling  /  rescaling  method  to  ‘fix’ 
the  spanwise  positions  of  longitudinal  vortices  (in  the  time  average)  that  are  generated  within  the 
upstream  boundary  layer.  Due  to  the  action  of  shearing  stresses,  these  structures  become  more 
elongated  near  the  wall,  and  because  of  the  absence  of  realistic  turbulence  in  the  logarithmic 
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Figure  23:  Surface  velocity  contours  (range:  -5  m/s  (dark)  to  35  m/s)  for  idealized  SBLI 

(LES/RANS  model) 

region  and  below,  the  ‘footprints’  of  the  outer- layer  structures  are  imposed  upon  the  near-surface 
velocity  field.  The  persistence  of  the  structures  for  this  particular  calculation  is  also  influenced 
by  the  relatively  coarse  mesh  (ranging  from  about  9  cells  per  boundary  layer  thickness  in  the 
wall-transverse  (X  and  Z)  directions  at  the  inflow  to  15  cells  per  boundary  layer  thickness  near 
the  exit). 

Because  of  the  dominant  effect  of  the  sidewall  separation  in  the  wind  tunnel,  it  is  difficult 
to  draw  any  direct  comparisons  between  the  idealized  and  wind-tunnel  cases.  It  can  be  noted, 
however,  that  the  area  occupied  by  low-momentum  fluid  is  similar  in  size  for  the  two  cases 
(Figures  22  and  14),  indicating  that  the  effect  of  the  crossflow  in  the  wind  tunnel  is  to  restrict  the 
growth  of  the  separation  region  in  the  spanwise  direction  while  enhancing  its  growth  in  the 
streamwise  direction. 

6.2.1 1.  Idealized  shock  /  boundary  layer  interaction  with  micro-VG  flow  control 


38 


Center  plane  Mach  number  contours  corresponding  to  the  idealized  oblique  shock  / 
boundary  layer  interaction  with  3  mm  micro  VG  control  are  shown  in  Figure  24  for  the  RANS 
(steady-state  contours)  and  LES/RANS  (instantaneous  contours)  models.  The  generation  of  an 
additional  oblique  shock  wave  by  the  micro  VG  is  of  note,  as  is  the  effect  of  the  incident  oblique 
shock  in  turning  the  vortex  pair  generated  by  the  micro  VG  more  toward  the  surface.  The 
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Figure  24:  Centerplane  Mach  number  contours  (idealized  interaction 

with  micro  VG  array) 

instantaneous  snapshot  indicates  a  more  energetic  eddy  structure  downstream  of  the  shock 
impingement  position.  This  may  result  from  turbulence  amplification  through  baroclinic  torque, 
bulk  compression,  shock  motion,  and  other  mechanisms. 

Surface  velocity  contours  are  shown  in  Figures  25  and  26  for  the  RANS  and  LES/RANS 
models,  respectively.  The  vortices  predicted  by  the  LES/RANS  model  transfer  more  high- 
momentum  fluid  to  the  surface  in  the  region  directly  behind  the  micro  VGs,  this  being  indicated 
by  lighter  contours.  As  a  consequence,  the  overall  separation  extent  is  reduced  for  the 
LES/RANS  calculation,  relative  to  the  RANS  calculation.  This  is  particularly  evident  for  the 
low-momentum  region  directly  behind  the  micro  VG  apex.  Comparing  Figures  25  and  26  with 
Figures  22  and  23,  it  is  seen  that  the  dominant  effect  of  the  micro  VG  array  is  to  break  up  the 
nominally  two-dimensional  separation  region  into  smaller  pockets.  The  micro  VG  array  appears 
more  effective  in  actually  eliminating  regions  of  low-momentum  fluid  in  this  idealized 
interaction  than  in  the  wind-tunnel  case.  The  vortices  generated  by  the  micro  VG  array  appear  to 
broaden  downstream  of  the  interaction  for  the  LES/RANS  model.  This  could  be  a  consequence 
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of  the  turbulence  amplification  mechanisms  noted  earlier,  and  the  effect  is  to  homogenize  the 
surface  fluid  in  the  spanwise  direction.  The  average  surface  velocity  at  the  exit  plane  is  about 


Figure  25:  Surface  velocity  contours  (range:  -5  m/s  (dark)  to  35  m/s)  for 
idealized  SBLI  with  micro  VG  control  (Menter  SST  RANS  model) 


30  m/s  for  the  LES/RANS  model  and  about  28  m/s  for  the  RANS  model,  but  the  spanwise 
distribution  of  velocity  for  the  LES/RANS  model  is  much  more  uniform.  The  region  of  low- 


Figure  26:  Surface  velocity  contours  (range:  -5  m/s  (dark)  to  35  m/s)  for 
idealized  SBLI  with  micro  VG  control  (LES/RANS  model) 


momentum  fluid  directly  behind  the  micro  VG  apex  is  smaller  in  the  LES/RANS  calculation, 
and  the  minimum  velocity  in  this  region  is  higher  (~25  m/s  versus  ~12  m/s).  The  fact  that  the 
vortices  generated  by  the  micro  VGs  continue  to  energize  the  boundary  layer  well  downstream  of 
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the  main  interaction  is  noteworthy,  as  it  indicates  the  potential  of  controlling  weaker  shock  / 
boundary  layer  interactions  generated  as  a  part  of  a  shock  train,  even  if  the  main  separation 
region  is  not  totally  eliminated. 

6.3.  Bleed  Control  of  Shock  /  Boundary  Layer  Interactions 

The  second  set  of  calculations  presented  in  this  work  utilizes  the  LES/RANS/IB 
methodology  to  simulate  the  effects  of  bleed-hole  arrays  in  controlling  shock-induced  flow 
separation.  Here,  an  entire  array  of  bleed  holes  (up  to  68),  along  with  portions  of  the  supporting 
flat  plate,  are  rendered  as  immersed  objects.  The  following  sections  are  excerpted  from  papers 
P-2  and  P-6  in  the  publication  list  (see  Section  9). 

6.3.1.  Experiment  details 

Experiments  involving  the  effects  of  boundary  layer  bleed  on  shock  interaction  control 
were  performed  at  NASA  Glenn  Research  Center  by  Willis,  et  al.  [26]  .The  experiments  were 
performed  at  a  nominal  Mach  number  of  2.46,  a  stagnation  temperature  of  292  K,  and  a 
stagnation  pressure  of  172.4  kPa.  In  the  configuration,  an  eight-degree  wedge  was  placed  so 
that  the  generated  oblique  shock  impinged  in  the  middle  of  a  9.52  cm  bleed  region  embedded 
within  the  wind-tunnel  lower  surface.  This  region  consisted  of  a  regular  array  of  one  hundred 
circular  holes  with  diameter  D  =  0.635  cm  and  length  L  =  0.635  cm,  oriented  90  degrees  with 
respect  to  the  plate  surface  The  bleed  flow  exited  into  a  plenum  chamber,  and  the  rate  of  mass 
flow  out  of  the  system  was  controlled  by  a  regulator  and  driven  by  a  vacuum  generated  by  a  450 
psi  air  ejector  system.  Both  the  plenum  pressure  and  overall  bleed  mass  flow  were  recorded  in 
the  experiments.  Data  collected  in  these  experiments  included  Pitot  pressure  surveys  at  different 
streamwise  locations,  wall  static  pressure  distributions,  and  Pitot  pressure  surveys  taken  in  the 
bleed-region  exit  plane,  all  as  a  function  of  the  bleed  mass  flow  rate.  In  the  experiments,  the 
bleed  mass  flow  rate  is  expressed  as  a  discharge  coefficient  Q,  representing  the  ratio  of  the  actual 
bleed  rate  to  the  ideal  value  obtained  under  choked-flow  conditions.  Data  was  obtained  for 
three  values  of  Q  -  0.00,  0.0342,  and  0.0685  -  as  well  as  for  a  baseline  case  where  the  bleed 
plate  was  replaced  by  a  solid  surface.  This  case  differs  from  the  Q  =  0.00  case,  as  the  latter 
allows  bleed  into,  and  injection  out  of  the  plenum.  Except  for  the  zero  bleed  case  (Q  =  0.0),  all 
other  cases  are  considered  in  this  study. 

In  addition  to  the  Willis,  et  al.  shock  /  boundary  layer  experiment  [26],  calculations  were 
also  performed  for  an  experiment  involving  Mach  2.45  boundary-layer  flow  over  a  bleed  plate 
[34],  This  experiment  focused  on  quantifying  the  discharge  coefficient  as  a  function  of  hole  size, 
hole  angle,  plenum  pressure,  and  free-stream  conditions.  The  hole  shape  is  the  same  as  in  [26] 
but  the  spacing  of  the  holes  within  the  array  is  different. 

6.3.2.  Computational  domain  and  calculation  details 

The  computational  domain  for  the  Willis,  et  al.  [26]  flat-plate  boundary-layer  simulations 
extends  from  X  =  -5.08  cm  to  X  =  17.78  cm  in  the  streamwise  direction,  from  Y  =  0  cm  to  Y  = 
7.62  cm  in  the  wall-normal  direction,  and  from  Z  =  -1.905  cm  to  Z  =  1.905  cm  in  the  spanwise 
direction.  The  bleed  plenum  for  this  case  ranges  from  X  =  -1.905  cm  to  X  =  9.525  cm  and  from 
Y  =  -7.62  cm  to  Y  =  -0.635  cm,  and  the  bleed-hole  region  extends  from  X  =  0  cm  to  X  =  7.62  cm 
and  from  Y  =  -0.635  cm  to  Y  =  0  cm.  To  provide  better  resolution  near  the  bleed  ports,  the 
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meshes  covering  the  bleed  region  are  refined  by  a  factor  of  two  in  the  Z  direction.  A  patched- 
mesh  boundary  condition  [35]  is  used  to  facilitate  the  coarse-to-fine  and  fme-to-coarse 
information  transfers  in  this  region.  The  X-Z  resolution  is  sufficient  to  cover  each  circular  bleed 
hole  with  a  10x10  square  mesh.  The  total  number  of  cells  is  3.24  M,  with  2.268  M  located 
above  the  flat  plate,  0.432  M  within  the  bleed-hole  region,  and  0.540  M  in  the  plenum.  The 
bleed  plate  itself  is  generated  as  a  combination  of  surfaces:  two  flat  plates  with  cut-out  holes  for 
the  top  and  bottom  portions  and  18  hollow  cylinders  (L  =  0.635  cm,  D  =  0.635  cm)  for  the  hole 
walls.  A  relatively  fine  IB  surface  mesh,  containing  1.32  M  points,  was  used  to  ensure  an 
accurate  distance-function  evaluation.  Figure  27  shows  the  immersed-surface  rendition  of  one 
circular  cylinder  (hole  wall)  and  the  surrounding  mesh. 

The  computational  domain  for  the 
Willis,  et  al.  [26]  shock  /  boundary-layer 
interaction  simulations  extends  from  X  =  - 
31.81  cm  to  X  =  25.0  cm  in  the  streamwise 
direction,  from  Y  =  0.0  cm  to  Y  =  21.5  cm 
in  the  wall-normal  direction  at  the  inflow 
plane,  and  from  Z  =  -5.08  cm  to  Z  =  5.08 
cm  in  the  spanwise  direction.  The  upper 
boundary  is  aligned  with  the  position  of  the 
8 -degree  shock  generator.  The  nominal 
mesh  spacing  in  the  X  and  Z  directions  is 
0.13  cm,  decreasing  to  half  of  this  value  in 
the  portions  of  the  mesh  covering  the  bleed 
region.  The  meshes  covering  the  bleed 
region  range  from  X  =  -2.23  cm  to  X  = 

13.01  cm  and  from  Y  =  -2.5  cm  to  Y  =  0.0 
cm.  The  plenum  chamber  extends  from  X  = 

-14.5  cm  to  X  =  26  cm  at  the  upper  end  (end 
of  bleed  region,  Y  =  -2.5  cm)  and  from  X  = 

-3.82  cm  to  X  =  15.82  cm  at  the  exit  (Y  =  - 
24.64  cm)  The  cross-section  of  the  plenum 
is  rectangular  till  about  Y  =  -9.88  cm  and 
then  tapered  to  the  dimensions  at  the  exit 
using  cubic  fits.  The  bleed  ports  for  this  case  are  the  same  size  as  in  the  flat-plate  study  and  the 
X-Z  resolution  is  sufficient  to  cover  each  port  with  10x10  square  cells.  The  total  number  of 
mesh  cells  is  21.504  M,  with  12.544  M  located  above  the  flat  plate,  1.792  M  located  in  the 
portion  of  the  domain  covering  the  bleed  holes,  another  1.792  M  in  the  region  between  the  bleed 
holes  and  the  plenum  and  5.376  M  located  in  the  plenum.  The  bleed-plate  surface  mesh 
contains  4.367  M  points  and  is  comprised  of  two  flat  sections  with  cut-out  holes  and  68 
individual  hollow  cylinders. 

No-slip,  adiabatic-wall  boundary  conditions  are  imposed  on  all  solid  surfaces  (excepting 
shock-generator  and  plenum  walls  where  a  slip  wall  condition  is  imposed)  for  both  cases.  A 
linear  power-law  is  used  to  define  the  near-surface  velocity  interpolation  for  the  flat-plate 
portions  of  the  immersed  bleed-plate  object,  as  the  mesh  spacing  near  the  wall  is  sufficient  to 
resolve  the  viscous  sublayer.  The  power-law  for  the  portions  representing  the  hole  walls  was 
varied  as  part  of  the  flat-plate  boundary-layer  study,  but  in  general,  a  fractional  power  law  (1/7  or 
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Figure  27:  Mesh  around  the  IB  rendition  of  a 
circular  bleed  hole 
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1/9)  is  employed,  as  the  mesh  resolution  near  the  hole  surfaces  does  not  extend  into  the  viscous 
sublayer.  A  subsonic  outflow  condition  which  imposes  the  experimental  mass  flow  rate  (using 
a  relaxation  technique)  is  used  at  the  plenum  exit.  Pressure  and  temperature  are  extrapolated  and 
density  is  calculated  based  on  equation  of  state.  The  normal  velocity  is  then  fixed  based  on  the 
value  of  mass  flux  imposed  and  the  density.  The  tangential  velocity  components  are  also  zeroed 
out.  Periodic  boundary  conditions  are  imposed  in  the  spanwise  directions.  It  should  be  noted  that 
the  experiments  of  Willis,  et  al.  [26]  used  aerodynamic  fences  to  isolate  the  interaction  region 
from  the  comer  flows  developing  in  the  wind  tunnel.  The  spanwise  extent  of  the  computational 
domain  does  not  extend  to  the  positions  of  the  fences  (+/-  7.94  cm),  and  as  such,  the  calculations 
represent  an  idealized  shock  /  boundary-layer  interaction  free  from  effects  of  mean  three- 
dimensionality.  For  the  flat-plate  boundary-layer  simulations,  flow  properties  are  extrapolated 
from  the  interior  on  the  upper  boundary  and  the  downstream  boundary.  For  the  shock  / 
boundary-layer  interactions,  slip-wall  conditions  are  applied  along  the  upper  boundary,  and 
properties  are  extrapolated  from  the  interior  at  the  downstream  boundary. 

In  all  cases,  a  flat-plate  simulation  at  the  test-section  conditions  was  performed  to 
determine  the  place  at  which  the  predicted  boundary  layer  properties  most  closely  matched  the 
experimental  data.  Flow  properties  at  the  inflow  plane  of  the  computational  domain  were  then 
extracted  from  the  flat  plate  solution  at  a  location  consistent  with  the  position  of  the  inflow  plane 
relative  to  the  measurement  location.  Convergence  of  the  RANS  calculations  was  ascertained 
using  three  measures:  relative  decrease  in  residual  norm,  constancy  of  surface  quantities,  and 
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Figure  28:  3-D  view  of  flow  over  a  perforated  plate  with  active 
suction  (Mach  number  contours  shown) 
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constancy  of  the  global  mass  flow  rate  and  bleed  mass  flow  rates.  The  bleed  mass  flow  rates 
were  evaluated  at  the  bleed-hole  entrance  and  exit  planes  as  well  as  the  plenum  exit  plane.  The 
LES/RANS  calculations  were  initiated  by  super-imposing  scaled  boundary-layer  fluctuations 
from  an  earlier  calculation  onto  part  of  a  converged  RANS  solution.  After  several  flow-through 
times  to  eliminate  initial  transients,  time-averaged  statistics  for  the  LES/RANS  calculations  were 
collected  over  a  minimum  of  seven  flow-through  times,  based  on  the  domain  length  and  free- 
stream  velocity. 

6.3.3.  Flat-plate  boundary  layer  simulations  with  bleed 

Simulations  of  Mach  2.5  flow  over  a  perforated  plate  were  performed  to  ascertain  the 
ability  of  the  methodology  to  predict  the  discharge  coefficient  as  a  function  of  plenum  pressure. 
The  Menter  BSL  RANS  model  was  used  for  these  calculations.  Mach  number  contours  in  Figure 
28  show  the  characteristic  pattern  of  flow  entrainment  into  the  bleed  holes.  The  immersed- 
boundary  object  that  represents  the  bleed  plate  is  also  shown  in  Figure  28.  A  large  pocket  of 
separated  flow  is  formed  at  the  upstream  edge  of  each  hole.  The  near-wall  fluid  accelerates  to 
supersonic  velocities  within  the  bleed  port  and  is  forced  to  change  direction,  moving  down  into 
the  plenum.  The  separation  region  forms  an  effective  converging  /  diverging  nozzle  that  first 
slows  the  supersonic  flow  down,  then  expands  it  to  higher  supersonic  velocities  before  entering 
the  plenum.  Figure  29  plots  the 
discharge  coefficient  (actual  mass  flow 
rate  divided  by  ideal  (isentropic)  mass 
flow  rate  at  sonic  conditions)  versus  the 
ratio  of  the  plenum  pressure  to  the 
upstream  stagnation  pressure.  The 
predictions  are  seen  to  be  in  reasonable 
agreement  with  the  experimental  data. 

The  predicted  mass  flow  rate  is  a  strong 
function  of  the  power  used  to  define  the 
band-cell  velocity  field  near  the  hole 
surfaces.  Using  a  value  of  k  =  1 
(consistent  with  a  linear  near-wall 
velocity  field  and  laminar  flow)  leads  to 
too  much  flow  separation  within  the 
bleed  holes  and  to  a  reduction  in  the 
discharge  coefficient.  Values  more 
consistent  with  a  turbulent  profile  (k  = 

1/7  or  k=l/9)  promote  more  flow 
attachment  and  provide  discharge 
coefficients  more  in  accord  with  the  experimental  data.  The  predicted  discharge  coefficients  at 
lower  plenum  pressures  are  somewhat  smaller  than  indicated  in  the  experiment.  The  effect  of 
reducing  the  X-Z  mesh  resolution  over  each  hole  from  10x10  cells  /  hole  to  5x5  cells  /  hole  is 
also  shown  for  the  lowest  plenum  pressure.  A  reduction  in  predicted  mass  flow  rate  of  about 
17%,  relative  to  the  baseline  10x10  resolution,  is  observed. 

6.3.4.  Oblique-shock  /  turbulent  boundary  layer  interaction  without  bleed 


Figure  29:  Flow  coefficient  versus  plenum 
pressure 
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A  snapshot  of  temperature  contours  along  the  X-Y  centerplane  is  shown  in  Figure  30  for 
the  baseline  shock  /  boundary  layer  interaction  without  bleed.  The  results  indicate  a  general 
thinning  of  the  boundary  layer  downstream  of  the  shock  reflection  as  well  as  local  ‘hot  spots’ 
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Figure  30:  Snapshot  of  centerplane  temperature  contours  for  obliquie  shock 

interaction  without  bleed 

(darker  contours)  caused  by  turbulent  separation  and  re-attachment  as  well  as  by  shock 
interaction.  The  shock  wave  is  strong  enough  to  induce  mild  separation  of  the  turbulent 
boundary  layer,  as  indicated  in  Figure  31,  a  plot  of  centerline  wall  pressure.  The  better 
agreement  with  experimental  data  is  provided  by  the  Menter  SST  RANS  model,  which  correctly 
predicts  the  degree  of  upstream  influence  shown  in  the  experiment.  The  LES/RANS  model  and 
Menter  BSL  RANS  model  (not  shown)  tend  to  under-predict  the  size  of  the  separation  region. 


Figure  31:  Wall  pressure  distributions:  shock  /  boundary  layer 
interaction  without  bleed 
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Centerline  Pitot-pressure  distributions  throughout  the  interaction  region  are  shown  in 
Figure  32  for  the  LES/RANS  and  Menter  SST  models.  The  first  trace  of  the  incident  oblique 
shock  is  captured  at  station  4,  with  the  sharp  pressure  change  at  ~  4  cm  from  the  wall.  Stations  5- 
8  (x  =  3.8,  29.2,  54.5,  80.0  mm)  form  the  core  of  the  interaction  of  the  oblique  shock  with  the 
boundary  layer.  This  is  expected  as  these  stations  are  within  the  bleed  plate  and  the  inviscid 
oblique  shock  impinges  ideally  at  the  center  of  the  bleed  region.  The  experimental  data  shows 
the  onset  of  separation  in  profile  5.  Good  agreement  is  generally  observed,  with  some 
discrepancies  in  the  wake-like  profile  noted  in  the  LES/RANS  predictions  at  the  X  =  3.8  mm  and 
X  =  29.2  mm  stations.  These  are  related  to  the  under-prediction  of  axial  separation  noted 
earlier.  In  the  recovery  region,  the  LES/RANS  and  RANS  predictions  generally  bracket  the 
experimental  data,  with  the  former  slightly  over-predicting  the  rate  of  recovery  of  the  boundary 
layer. 
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Figure  32:  Pitot  pressure  distributions:  shock  /  boundary  layer  interaction  without 

bleed 


6.3.5.  Oblique-shock  /  turbulent  boundary  layer  interaction  with  bleed 

The  cases  detailed  in  this  section  involve  suction  through  a  68-hole  perforated  plate, 
shown  in  Figure  33  as  the  zero  iso-surface  of  the  signed  distance  function.  Two  different  bleed 
rates,  referred  to  henceforth  as  ‘full  bleed’  (maximum  bleed)  and  ‘half  bleed’,  which  is  basically 
half  of  the  maximum  bleed,  were  simulated.  Both  RANS  and  LES/RANS  simulations  were 
performed  for  the  full  bleed  case  while  only  the  LES/RANS  model  was  used  for  the  half  bleed 
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simulation.  Both  the  simulations  were  initialized  with  plenum  pressures  which  correspond  to 
those  in  the  experiment  for  the  respective  bleed  rates  (5534  Pa  for  full  bleed  and  17413  Pa  for 
half  bleed).  In  the  half-bleed  case,  the  bleed  mass  flow  rates  are  much  more  unsteady  due  to 
periodic  blowing  and  suction  through  the  holes  located  upstream  of  the  shock  impingement 
position. 


a)  half  bleed 


b)  full  bleed 

Figure  33:  Snapshot  of  temperature  contours  with  inset  bleed  plate:  shock  / 
boundary  layer  interaction  with  bleed 
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Temperature  contours  shown  in  Figure  33  reveal  that  the  bleed  array  produces  localized 
wave  patterns  associated  with  the  formation  of  expansion  waves  as  the  near-wall  fluid  moves 
into  each  bleed  port  and  with  the  formation  of  a  re-compression  shock  located  at  the  downstream 
edge  of  the  bleed  port.  This  response  occurs  for  all  the  bleed  ports  in  the  full  bleed  case  but  only 
in  the  most  downstream  (last  two)  ports  for  the  half  bleed  case.  A  general  thinning  of  the 


a)  half  bleed 


b)  full  bleed 

Figure  34:  Wall  pressure  distributions:  shock  /  boundary  layer 
interaction  with  bleed;  (*)  -  averaged  at  pressure  tap  locations 
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boundary  layer  within  and  downstream  of  the  bleed  region  is  evident,  as  is  a  re-compression 
shock  located  just  downstream  of  the  end  of  the  bleed  region.  Overall  the  effects  are  more 
pronounced  for  the  full  bleed  case.  Higher  temperatures  are  noted  near  the  vicinity  of  the  shock 
impingement  position  for  the  half  bleed  case.  These  are  associated  with  a  recirculation  region 
formed  due  to  the  shock  impingement  and  enlarged  somewhat  by  air  blowing  out  from  the 
plenum.  The  free-jet  expansion  of  the  bleed-hole  fluid  as  it  enters  the  plenum  is  less  pronounced 
for  in  the  half-bleed  case,  due  to  the  higher  plenum  pressure. 

Centerline  surface  pressure  distributions  in  Figure  34  provide  more  evidence  of  the 
localized  expansion  /  compression  events  induced  by  the  suction  effect.  The  pressure  levels  at 
the  surface  plane  over  each  bleed  port  are  included  in  the  computational  distributions  to  provide 
an  indication  of  the  local  flow  patterns.  For  the  full  bleed  simulation,  the  LES/RANS  and  RANS 
predictions  are  very  similar,  and  when  averaged  to  the  experimental  data  locations,  are  in 
reasonable  agreement  with  the  measurements.  Surface  pressure  predictions  for  the  half  bleed 
case  also  compare  reasonably  well  with  the  data.  It  is  interesting  to  note  that  the  localized 
compression/expansion  effects  are  relatively  absent  in  the  upper  part  of  the  interaction  region  for 
the  half  bleed  case  compared  to  the  maximum  bleed  scenario.  This  is  due  to  the  absence  of  active 
suction  through  the  bleed  ports  in  this  region  for  the  lower  bleed  rate.  In  Ref.  [26],  some 
discussion  is  made  of  the  fact  that  the  measured  wall  pressure  for  the  maximum  bleed  simulation 
is  higher  in  some  places  than  the  value  obtained  through  oblique  shock  theory.  The  authors 
argue  that  situations  might  arise  where,  due  to  the  suction  effect,  flow  is  directed  into  the  static 
tap,  causing  it  to  read  a  higher  value.  Another  explanation,  furthered  by  the  computational 
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Figure  35:  Pitot  pressure  distributions:  shock  /  boundary  layer  interaction  with  half¬ 
bleed 
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results,  might  be  that  some  of  the  static  taps  are  very  close  to  a  bleed  hole  and  as  such,  the  static 
taps  measure  an  average  value  associated  with  the  rapid  compression  /  expansion  process 
occurring  at  the  downstream  edge  of  each  bleed  hole.  It  is  clear  that  the  bleed  effects  are  highly 
localized,  at  least  for  this  particular  array. 

Pitot-pressure  distributions  throughout  the  interaction  region  are  shown  in  Figures  35  and 
36  for  the  half-and  full-bleed  simulations.  The  model  predictions  are  very  similar  and  agree 
well  with  the  experimental  data.  The  LES/RANS  model  produces  a  slightly  fuller  Pitot-pressure 
profile  near  the  wall  downstream  of  the  bleed  region  for  both  the  bleed  rates.  This  result  is 
consistent  with  that  evidenced  for  the  case  without  bleed,  in  that  the  LES/RANS  model  generally 
provides  a  more  rapid  recovery  of  the  boundary  layer  downstream  of  a  shock-induced 


*  Expt.  -  RANS  -  LES/RANS 

x  =  80.0  mm  x=  115  mm  x=  130  mm  x=  156  mm  x=  182  mm  x=  207  mm  x  =  232  mm 


0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8 


x=  -98.0  mm  x  =  -82.0  mm  x  =  -47.0  mm  x=  -22.0  mm  x=  3.8  mm  x=  29.2  mm  x=  54.6  mm 


0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8  0.2  0.4  0.6  0.8 


Pt/PO 

Figure  36:  Pitot  pressure  distributions:  shock  /  boundary  layer  interaction  with  full 

bleed 

disturbance.  Predictions  within  in  the  interaction  region,  specifically  at  the  80  mm  station,  are 
not  quite  as  good,  especially  for  the  half  bleed  case.  This  is  possibly  due  to  the  difference  in  the 
structure  of  the  separation  shock  between  experiment  and  that  predicted  by  computation. 

Pitot-pressure  distributions  extracted  at  the  bleed-hole  exit  plane  (Figure  37)  for  the  full 
bleed  case  indicate  that  the  computations  under-predict  the  measured  peak  pressure  level  at  all 
axial  stations,  with  the  LES/RANS  model  providing  slightly  better  agreement.  This  trend  is 
consistent  with  the  under-prediction  of  the  overall  bleed  mass  flow  rate  mentioned  above.  It  is 
possible  that  these  results  would  improve  with  additional  X-Z  mesh  refinement  over  the  holes. 

The  bleed  mass  flux  distributions  over  the  entire  array  for  the  full-  and  half-bleed 
LES/RANS  simulations  are  displayed  in  Figure  38.  Of  particular  note  is  the  increase  in  bleed 
mass  flux  in  the  region  downstream  of  the  shock-impingement  location.  This  is  the  expected 
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trend,  as  the  pressure  level  is  higher  in  this 
region.  For  the  full-bleed  case,  the 
suction  effect  is  much  less  for  bleed  ports 
located  upstream  of  the  shock 
impingement  location.  The  mass  fluxes  in 
this  region  are  similar  to  those  found  in 
the  flat-plate  boundary  layer  simulations 
discussed  above.  Some  localized 
entrainment  of  bleed-hole  fluid  into  the 
boundary  layer  is  noted  for  the  full-bleed 
case,  but  a  much  stronger  blowing  effect  is 
present  for  the  half-bleed  case.  This  is 
due  to  a  pressure  differential  within  the 
plenum  caused  by  the  preferential 
injection  of  mass  downstream  of  the  shock 
impingement  location.  Mach  number 
contours  and  streamlines  in  Figure  39 
indicate  that  the  initial  entrainment  of 
fluid  into  the  bleed  port  occurs  through  a 
supersonic  expansion  that  is  terminated  by 
a  ‘barrier  shock’  [13]  oriented  normal  to  the  flow  direction.  The  lowermost  portion  of  this  shock 
wave  interacts  with  the  fluid  near  the  upstream  edge  of  the  bleed  port,  creating  a  region  of 


a)  half  bleed  b)  full  bleed 

Figure  38:  Bleed  mass  flux  distributions 


separation  that  acts  to  reduce  the  cross-section  area  occupied  by  the  core  fluid.  The  growth  of 
this  region  creates  an  area  minimum,  which  allows  the  subsonic  flow  behind  the  barrier  shock  to 
accelerate  to  supersonic  speeds  at  the  bleed-hole  exit.  This  under-expanded  jet  then  relaxes  to 
the  plenum  pressure  level  through  a  series  of  expansion  and  compression  waves.  The  uppermost 
portion  of  the  barrier  shock  propagates  out  of  the  bleed  port.  The  fluid  that  passes  through  this 
part  of  the  barrier  shock  is  compressed  to  levels  in  excess  of  the  average  wall  pressure  and  is 


Figure  37:  Pitot  pressure  profiles  at  exit  of  bleed 
holes  with  full  bleed 


51 


0.005 


Barrier  shock 


0 


-0.005 

I 

*-0.01 


-0.015 


-0.02 


Sonic  line 


2.40773 

2.1141 

1.82048 

1 .52685 

1.23323 

0.9396 

0.645975 

0.35235 

0.058725 


0.07 


0.08 

x[m] 


J _ I _ L 

0.09 


i 


Figure  39:  Mach  number  contours  in  bleed  hole 


then  expanded  as  it  turns  toward  the  flat-plate  surface.  This  provides  the  mechanism  for  the 
large  spikes  in  the  wall  pressure  distributions  shown  in  Figure  34. 


6.3.6.  Reynolds-stress  evolution  and  turbulence  structure 


Turbulence  amplification  is  a  well-known  characteristic  of  shock  /  boundary  layer 
interactions.  Bulk  compression,  streamline  curvature,  shock  oscillation,  and  the  dynamics  of 
separation  /  re-attachment  all  contribute  to  the  amplification  effect.  [36]  Large  bleed  rates  as 
used  in  the  Willis,  et  al.  experiments  can  essentially  remove  turbulent  separation  but  also  induce 
other  disturbances  that  could  potentially  amplify  turbulence  Figure  40  plots  span-averaged 
contours  of  resolved  turbulence  kinetic  energy  and  resolved  Reynolds  shear  stress  ( pu'v' )  for  the 
LES/RANS  cases  with  and  without  bleed.  For  the  case  without  bleed,  large  Reynolds-stress 
values  are  found  near  the  time-averaged  position  of  the  separation  shock  and  downstream  of  the 
re-attachment  location.  Further  downstream,  the  stresses  begin  to  relax  to  levels  similar  to,  but 
larger  than  those  in  the  incoming  boundary  layer.  With  full  bleed,  the  separation  /  re-attachment 
response  is  eliminated,  and  very  little  Reynolds-stress  amplification  is  observed.  The  half-bleed 
results  indicate  that  some  initial  amplification  takes  place  near  and  upstream  of  the  shock- 
impingement  position,  as  there  is  no  significant  suction  effect  in  this  region.  Further 
downstream,  the  suction  effect  suppresses  the  shock  system  dynamics,  and  the  Reynolds  stresses 
again  begin  to  relax  to  values  similar  to  those  in  the  upstream  boundary  layer.  The  wave  patterns 
generated  by  flow  into  the  bleed  holes  generate  local  disturbances  in  the  Reynolds  stresses,  but 
these  propagate  out  of  the  boundary  layer.  This  result  indicates  that  the  dominant  source  of 
Reynolds-stress  amplification  in  impinging  shock  /  boundary-layer  interactions  may  be  the 
dynamics  of  the  separation  region,  which  triggers  the  motion  of  the  separation  and  re-attachment 
shocks. 

Figure  41  shows  iso-surfaces  of  the  swirl  strength  (10000  s'1)  for  the  interactions  with 
and  without  full  bleed.  Visually,  the  interactions  are  similar,  displaying  a  general  increase  in 
structural  content  and  the  presence  of  more  smaller-scale  structures  within  and  downstream  of 
the  interaction  region.  A  more  quantitative  measure,  shown  in  Figure  42,  is  obtained  by 


52 


calculating  the  probability  density  function  of  the  logarithm  of  the  swirl  strength  at  different 
streamwise  locations.  The  most  probable  value  is  indicated  as  the  peak  in  the  distribution,  and  a 
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Figure  40:  Evolution  of  resolved  turbulent  kinetic  energy  (left)  and  resolved  Reynolds 
shear  stress  (right)  through  the  interaction  region 
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Figure  41:  Iso-surfaces  of  swirl  strength  (10000  s-1)  for 
interactions  with  full  bleed  and  no  bleed 

shift  in  the  distribution  to  the  right  implies  an  increase  in  the  population  of  more  tightly-wrapped 
(fmer-scale)  vortical  structures.  This  can  also  imply  a  general  decrease  in  turbulence  length 
scales.  For  the  interaction  without  bleed,  the  most  probable  swirl  strength  value  increases  by 
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probability  density 


about  an  order  of  magnitude  as  the  flow  separates,  then  re-attaches.  In  the  recovery  region,  the 
most  probable  swirl-strength  value  decreases,  and  the  distribution  approaches  that  of  the 
upstream  boundary  layer.  With  bleed,  the  distributions  do  not  deviate  as  much  in  the  interaction 
region,  though  the  trend  of  an  increasing  population  of  finer-scale  structures  still  holds.  The 
swirl-strength  distributions  in  the  recovery  region  are  close  to  those  in  the  incoming  boundary 
layer,  again  showing  that  the  bleed  effect  accelerates  the  recovery  of  the  boundary  layer  to  a  new 
equilibrium  state  following  the  shock  reflection. 


Figure  42:  Evolution  of  swirl-strength  probability  density  distributions 


6.4.  Mesoflap  Control  of  Shock  /  Boundary  Layer  Interactions 

The  third  set  of  calculations  presented  in  this  work  applies  the  LES/RANS/IB 
methodology  to  examine  shock  /  boundary  layer  interaction  control  using  arrays  of  mesoflaps, 
which  are  rendered  as  immersed  bodies  and  can  deform  dynamically  in  response  to  time-varying 
pressure  loads.  We  consider  quasi-steady  simulations  using  a  2-D  RANS  solver,  similar  to 
Wood,  et  al.  [32]  as  well  as  fully  dynamic  simulations  using  the  LES/RANS  model. 

6.4.1.  Experiment  details 

Experiments  involving  mesoflap  arrays  having  different  thicknesses,  layouts  and 
arrangement  to  control  SBLIs,  have  been  performed  at  the  University  of  Illinois  [6-8].  The 
experiments  were  performed  at  a  nominal  Mach  number  of  2.45,  a  stagnation  temperature  of  300 
K,  and  a  Reynolds  number  /  meter  of  57xl06.  The  test  section  in  the  wind  tunnel  is  50.8  mm 
high  and  50.8  mm  wide.  Data  collected  in  these  experiments  includes  centerline  axial  velocity 
data  obtained  at  various  streamwise  stations  using  one-component  laser  Doppler  velocimetry 
(LDY)  and  wall  static  pressure  distributions.  The  velocity  data  has  been  used  for  generating 
profiles  of  streamwise  normal  stress  intensity  and  boundary  layer  integral  properties.  This  work 
focuses  on  two  of  the  UIUC  experiments: 
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•  An  oblique  shock  /  turbulent  boundary  layer  interaction  at  Mach  2.45  induced  by  a  8- 

degree  wedge  placed  on  the  upper  wall  of  the  wind  tunnel  (LDV  velocity  data, 
surface  pressure  measurements) 

•  The  same  oblique  shock  /  turbulent  boundary  layer  interaction  with  the  addition  of  an 

array  of  six  third-generation  [7]  mesoflaps  of  thickness  =  76.2  pm  in  the  shock 
interaction  region  (LDV  velocity  data,  surface  pressure  measurements). 

The  mesoflap  system  consisted  of  a  steel  frame  with  supporting  spars  with  the  mesoflap 
array  made  of  nickel-titanium  (nitinol)  epoxied  to  it.  The  mesoflaps  were  4.63  mm  long,  spanned 
about  %  of  the  wind  tunnel  and  had  a  thickness  of  76.2  pm.  The  mesoflap  system  enclosed  a 
cavity  which  was  44.5  mm  wide  and  19.1  mm  deep.  The  leading  edge  of  the  first  mesoflap 
(counted  from  upstream  to  downstream)  was  located  at  -22.25  mm  from  the  center  of  the  cavity 
(which  was  also  the  shock  impingement  location).  This  location  has  been  denoted  as  x0  in  the 
experiment  and  also  in  this  work.  For  convenience,  the  computational  domain  used  in  this  work 
(detailed  in  the  next  section)  sets  x0  =  0.  The  leading  edge  of  the  shock  generator  was  placed  at 
19.45q  upstream  of  xq,  which  was  also  the  location  where  the  inflow  boundary  layer  properties 
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Figure  43:  3-D  views  of  mesoflap;  a)  deflected  mesoflap  array  with  flow  direction  indicated,  b)  a 
single  mesoflap  (entire  width  in  z  direction  not  shown  for  compactness) 

were  calculated  in  the  experiment.  In  the  computational  domain,  the  shock  generator  position 
was  adjusted  such  that  the  inviscid  oblique  shock  impingement  location  was  approximately  at  xo. 
The  boundary  layer  thickness  at  the  inflow  plane,  So  was  determined  as  4.0  mm  by  fitting  the 
velocity  data  to  the  experimental  curve  fit  due  to  Sun  and  Childs  [37],  The  compressible 
displacement  thickness  and  momentum  thickness  were  also  calculated  from  the  curve  fit  as 
0.275o  and  0.0755o  respectively.  Velocity  profiles  were  collected  at  streamwise  stations  of  x  =  - 
19.48omm,  x  =  -0.12  So,  x  =  2.26  8o,  x  =  8.738o,  x  =  15.08o  for  the  SBLI  without  any  control.  For 
the  case  with  the  mesoflap  array  in  place,  velocity  data  was  collected  at  locations  downstream  of 
the  mesoflap  array  at  x  =  8.738o,  x  =  13.48o  [7]. 

6.4.2.  Computational  domain  and  calculation  details 
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The  primary  domain  (that  containing  the  shock  generator  and  incoming  boundary  layer) 
extends  from  X  =  -99.05  mm  to  X  =  60.65  mm  in  the  streamwise  direction  and  from  Y  =  0  to  Y 
=  58.0  mm  in  the  wall  normal  direction  at  the  start  of  the  domain,  reducing  to  ~  39.8  mm  at  the 
trailing  edge  of  the  shock  generator.  For  the  2-D  quasi-steady  FSI  simulations,  the  spanwise 
width  of  the  domain  is  one  mesh  cell.  For  the  LES/RANS  simulations,  the  computational  domain 
extends  from  Z  =-19.4  mm  to  Z  =19.4  mm  in  the  spanwise  direction.  The  streamwise  and 
spanwise  mesh  spacing  is  ~  5©/15,  and  the  minimum  mesh  spacing  in  the  wall-normal  direction  is 
0.005  mm.  For  the  cases  with  the  mesoflap  array  in  place,  another  domain  for  the  plenum  was 
added.  The  plenum,  which  is  covered  on  its  top  by  the  mesoflap  array,  extends  from  X  =  -22.25 
mm  to  X  =  22.25  mm  in  the  streamwise  direction  and  from  Y  =  0.  mm  to  Y  =  -19.1  mm  in  the 
wall  normal  direction.  The  spanwise  extent  of  the  cavity  is  same  as  for  the  corresponding 
primary  domain.  The  number  of  cells  used  for  the  upper  domain  in  the  LES/RANS  computations 
is  ~18.66M  with  the  plenum  containing  an  additional  ~1.55M  cells.  The  array  of  mesoflaps 
(including  the  spars)  is  rendered  as  a  collection  of  IB  surfaces.  Figure  43  shows  a  single  flap 
structure  with  the  fixed  spar.  The  number  of  points  required  in  rendering  the  array  structure  is 
791453. 

No-slip,  adiabatic-wall  boundary  conditions  are  imposed  on  all  solid  surfaces  (excepting 
the  upper  wall  including  the  shock-generator  and  plenum  walls  where  a  slip  wall  condition  is 
imposed)  for  both  cases.  The  value  of  power  law  chosen  in  different  parts  of  the  mesoflap  is  as 
follows  -  k  =  1  for  the  flap  upper  and  lower  surface  and  k  =  1/9  for  all  other  surfaces.  A  linear 
power-law  is  used  to  define  the  near-surface  velocity  interpolation  for  surfaces  in  which  the 
mesh  spacing  (normal)  near  the  surface  is  sufficient  to  resolve  the  viscous  sub  layer.  Periodic 
conditions  are  used  in  the  spanwise  directions.  For  the  shock  /  boundary-layer  interactions,  slip- 
wall  conditions  are  applied  along  the  upper  boundary,  and  properties  are  extrapolated  from  the 
interior  at  the  downstream  boundary.  Boundary  conditions  for  the  flat  plate  boundary  layer  study 
were  same  (on  the  lower  wall,  upper  wall  and  downstream  boundary)  as  for  the  SBLI  studies. 

To  determine  the  inflow  conditions  for  all  the  SBLI  studies,  a  flat-plate  simulation  at  the 
test-section  conditions  was  first  performed  to  determine  the  place  at  which  the  predicted 
boundary  layer  properties  most  closely  matched  the  experimental  data  at  X  =  -19.4  do.  Flow 
properties  at  the  inflow  plane  of  the  computational  domain  were  then  extracted  from  the  flat 
plate  solution  at  a  location  consistent  with  the  position  of  the  inflow  plane  relative  to  the 
measurement  location.  Convergence  of  the  RANS  calculations  was  ascertained  using  three 
measures:  relative  decrease  in  residual  norm,  constancy  of  surface  quantities,  and  constancy  of 
the  global  mass  flow  rate.  For  the  quasi-steady  FSI  simulation,  the  convergence  of  the  mesoflap 
position  was  monitored  using  changes  in  the  computed  displacements  of  the  flaps.  The 
LES/RANS  calculations  were  initiated  by  super-imposing  scaled  boundary-layer  fluctuations 
from  an  earlier  calculation  onto  part  of  a  converged  RANS  solution.  After  several  flow-through 
times  to  eliminate  initial  transients,  time-averaged  statistics  for  the  LES/RANS  calculations  were 
collected  over  a  minimum  of  seven  flow-through  times,  based  on  the  domain  length  and  free- 
stream  velocity.  The  LES/RANS  calculation  involving  the  mesoflap  array  was  run  for  a  longer 
duration  ~12  flow  through  times.  This  was  done  to  discern  a  pattern  in  the  mesoflap 
displacements  and  to  determine  the  dominant  response  frequencies. 
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6.4.3.  Oblique  shock  /  boundary  layer  interaction  without  mesoflap  array 


A  3-D  view  showing  snapshots  of  pressure  contours  on  a  streamwise  X-Y  plane  and  axial 
velocity  on  a  near  surface  plane  (X-Z)  is  shown  in  Figure  44  for  the  LES/RANS  model.  The  iso¬ 
surface  of  zero  axial  velocity,  indicating  the  separation  bubble,  is  colored  in  blue.  The  figure 


Figure  44:  Snapshot  of  contour  plots  of  pressure  on  an  X-Z  plane 
(close  to  surface)  and  axial  velocity  an  X-Y  plane;  iso-surface  of 
zero  axial  velocity  shown  in  blue  on  the  near  surface  plane 


shows  the  impinging  and  reflected  shock  structure  showing  that  the  zero  axial  velocity  iso¬ 
surface  originates  approximately  at  the  location  of  separation  shock.  The  separation  bubble  has 
‘bumps’  and  ‘troughs’  which  can  be  discerned  when  seen  closely.  These  are  due  to  streaks  of 
high  and  low-momentum  fluid  in  the  incoming  boundary  layer  close  to  the  surface.  Figure  45 
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Figure  45:  Snapshot  of  Mach  number  contours  showing  in 
white  the  zero  axial  velocity  boundary 

shows  a  2-D  close  up  view  of  the  separation  region  with  a  white  line  indicating  the  location  of 
zero  axial  velocity. 
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Figure  46  shows  the  centerline  wall  pressure  distribution  compared  with  experimental 
data  for  the  oblique  shock  interaction  with  no  control.  The  Menter  SST  RANS  model  appears  to 
over-estimate  the  level  of  upstream  influence  of  the  interaction  (indicated  by  the  separation- 
shock  pressure  rise).  The  LES/RANS  model  provides  slightly  better  predictions.  In  the 
experiment,  a  gradual  rise  in  pressure  is 
observed  upstream  of  the  interaction 
region  from  about  ({X - X0 ) / S(]  =-10, 

which  has  been  attributed  to  effects  of 
reduced  free  stream  Mach  number  due 
to  boundary  layer  growth  in  the  test 
section.  This  is  not  reflected  in  the 
numerical  data  and  might  be  due  to  the 
fact  the  boundary  layer  on  the  upper 
wall  (treated  as  a  symmetry  boundary  in 
computations)  is  not  resolved.  The 
computations  over-predict  the  pressure 
level  downstream  of  the  impingement 
position.  It  is  not  clear  why  this  might  ’  ’  ’  (x-x„)/80 

be  happening.  Beyond  this  region  the  Figure  46:  Centerline  wall  pressure  for  SBLI 

expansion  fan  generated  at  the  trailing  without  control 

edge  of  the  shock  generator  hits  the 
lower  wall,  which  is  reflected  by  the 

drop  in  the  surface  pressure.  The  predictions  in  this  region  are  accurate  for  both  computations. 

Figure  47  shows  experimental  and  computed  centerline  axial  velocity  profdes  at  different 
streamwise  locations.  The  numerical  predictions  of  axial  velocity  at  the  inflow 
((X-X0)/S0  =  -19.4)  are  in  good  agreement  with  experimental  data,  with  both  RANS  and 

□  Expt.  -  RANS  - LES/RANS 


Figure  47:  Centerline  axial  velocities  for  shock/boundary  layer  interaction  without  control 
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LES/RANS  profiles  being  somewhat  more  full  in  the  near-wall  region  (y  1 80<QA).  In  the 
interaction  region  (  (X-X0)  /  S0  =  -0.12,2.26  )  the  numerical  predictions  show  significant 

deviations  from  experimental  data.  The  wake-like  profile  present  in  the  experiment  is  less 
evident  in  the  computed  profiles,  more  so  for  the  LES/RANS  predictions.  This  is  most  evident  at 
the  (X  -  X0 )  /  S0  =  2.26  station.  This  is  suggestive  of  differences  in  the  shape  and/or  location  of 

the  separation  bubble  between  the  experiment  and  the  computations.  The  experimental  profile  is 
suggestive  of  a  thicker  separation  bubble.  A  possible  reason  for  this  might  be  the  fuller  near-wall 
velocity  profile  predicted  at  the  inflow,  which  would  be  more  resistant  to  separation.  The  last 
two  stations  are  downstream  of  the  interaction,  and  the  computational  predictions  are  closer  to 
the  LDV  data.  Here,  the  LES/RANS  method  tends  to  over-predict  the  recovery  rate  of  the  inner 
part  boundary  layer,  compared  to  the  RANS  model.  This  trend  has  been  observed  in  some 
previous  applications  of  the  LES/RANS  method  in  compression-ramp  interactions  [15]  and  may 
be  indicative  of  higher  turbulence  intensity  near  the  wall. 

Figure  48  compares  predicted  streamwise  turbulence  intensity  profiles  (resolved  stresses 
for  the  LES/RANS  model)  at  different  streamwise  locations  with  experiment.  Three  sub-plots  are 
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Figure  48:  Streamwise  turbulence  intensity  profiles  at  various  streamwise 
locations;  a)  inflow,  b)  interaction  region,  and  c)  downstream 

developed  by  categorizing  the  streamwise  stations  into  inflow  (a),  interaction  (b)  and 
downstream  (post-interaction;  c)  zones.  The  turbulence  intensity  levels  are  on  the  average  under¬ 
predicted.  At  the  inflow  station,  the  free-stream  turbulence  level  is  rather  high  ~2%  compared  to 
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the  calculation.  In  the  interaction  region,  the  location  of  maximum  turbulence  intensity  shifts 
away  from  the  wall.  Within  the  separation  region  ((X  -  X0)/ <5 ()  =  -0.12),  the  location  of  the 

peak  is  under-predicted  but  is  in  better  agreement  with  experiment  at  stations  further 
downstream.  The  turbulence  intensity  level  is  under-predicted  at  the  station  just  downstream  of 
the  re-attachment  position  ( (X  -  X0)/  S0  =2.26).  The  turbulence  intensity  profiles  downstream 
of  the  interaction  show  better  agreement  with  experimental  data,  with  excellent  agreement 
observed  at  the  (X - X0) / S0  =  8.73  station.  .  Overall  the  LES/RANS  predictions  fare  reasonably 

well  compared  with  experiment,  with  better  results  obtained  upstream  and  downstream  of  the 
interaction.  This  is  again  suggestive  of  differences  in  separation-bubble  geometry  and  possibly 
interaction  dynamics  in  that  region. 


Figure  49:  Mach  number  contours  for  the  converged  configuration  of  the  mesoflap  array 


6.4.4.  Flow  structure  of  oblique  shock  /  boundary  layer  interaction  with  mesoflap  array  -  2-D 

quasi-steady  FSI  modeling  (RANS) 

Two  different  approaches  have  been  used  to  account  for  the  effects  of  meso-flap  control: 
a  2-D  simulation  which  uses  a  quasi-steady  modeling  for  the  structure  and  uses  a  RANS 
turbulence  (Menter  SST)  closure,  and  a  3-D  simulation  which  employs  a  dynamic  modeling  for 
the  structure  and  uses  the  hybrid  LES/RANS  turbulence  closure. 

Figure  49  shows  the  Mach  number  contours  for  the  converged  configuration  of  the 
mesoflaps  as  modeled  using  the  quasi-steady,  2-D  approach.  Fluid  is  entrained  into  the  cavity 
through  the  slots  furthest  downstream  of  the  shock  impingement  point  and  is  injected  from  the 
cavity  into  the  boundary  layer  through  the  slots  furthest  upstream  of  the  shock  impingement 
point.  The  mesoflaps  are  numbered  from  the  upstream  end,  so  mesoflap  1  is  the  farthest 
upstream  and  mesoflap  6  is  farthest  downstream.  As  can  be  observed,  mesoflaps  1  and  2  and  5 
and  6  are  primarily  responsible  for  the  flow  into  and  out  of  the  cavity.  Mesoflaps  3  and  4  deflect 
very  little,  allowing  only  minimal  fluid  exchange  between  the  upper  domain  and  the  plenum. 
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This  trend  is  very  consistent  with  that  observed  in  the  experiments  [6,7],  Figure  50  shows  the 
maximum  deflection  (which  is  also  the  edge  deflection)  of  each  mesoflap  at  each  structural 
iteration.  To  improve  the  code’s  response  in  the  initial  transient  periods,  a  first  order  spatial 
discretization  was  used  for  fluid  solver  to  begin  with,  and  when  the  deflection  pattern  showed  a 
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Figure  50:  History  of  maximum  deflection  of  each  mesoflap  with  iteration 

reasonable  convergence,  the  higher-order  PPM  discretization  was  turned  on.  There  is  a  sharp 
change  in  the  deflection  pattern  of  the  upstream  flaps  at  around  the  34th  structural  iteration  which 
roughly  corresponds  to  this  change  in  the  order  of  accuracy.  It  is  interesting  to  note  that  the 
upstream  mesoflaps  (1-2)  tend  to  converge  faster  than  the  downstream  (5-6)  ones. 


Figure  51:  Mach  number  contours  for  shock/boundary  layer  interaction 
with  control  by  mesoflap  array 
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6.4.5.  Flow  structure  of  oblique  shock  /  boundary  layer  interaction  with  mesoflap  array  -  3-D 

dynamic  FSI  modeling  (LES/RANS) 


Figure  51  shows  a  3-D  view  of  the  deflected  mesoflap  array  with  a  snapshot  of  Mach 
number  contours  shown.  The  trace  of  the  leading  edge  shock  can  be  seen  on  the  X-Y  plane,  with 
a  low  Mach  number  region  being  present  over  a  large  portion  of  the  mesoflap  array.  A  sequence 
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Figure  52:  Snapshots  of  Mach  numbers  contours;  frames  are  spaced  at 
intervals  of  600  iterations  ~  0.09  ms 

of  Mach  number  snapshots  is  shown  in  Figure  52  to  illustrate  the  dynamics  of  the  flow  with  the 
mesoflap  array  in  place.  The  sequence  spans  ~4200  iterations  which  covers  about  0.63  ms  in 
physical  time.  It  can  can  be  seen  that  injection  and  bleed  happen  simultaneously  at  the  upstream 
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and  downstream  mesoflaps,  respectively.  The  angle  at  which  the  fluid  is  injected  from  the  cavity 
into  the  core  flow  or  bleeds  into  the  cavity  from  the  core  flow  is  dependant  on  flap  configuration. 
As  observed  in  the  quasi-steady  simulation,  it  is  evident  that  the  intermediate  flaps  (3,  4)  do  not 
play  a  major  role  in  the  transpiration  process.  Angled  injection  at  the  most  upstream  flap  seems 
to  energize  the  flow  (leading  to  a  reduction  in  the  low  Mach  number  region)  whereas  at  the 
second  flap,  injection  tends  to  swell  up  the  separation  bubble  forming  just  downstream  of  it.  The 
maximum  amount  of  bleed  occurs  at  the  most  downstream  mesoflap. 


Figure  53:  Wall  static  pressure  history  with  the  time 
averaged  value  shown  in  red 

Figure  53  shows  the  history  of  wall  pressure  spanning  4200  iterations  which  is  about  0.63 
ms,  with  data  spaced  at  200  iterations.  The  figure  also  shows  the  location  of  the  mesoflap  array 
in  the  lower  part  of  the  frame  to  help  identify  patterns  in  the  time-averaged  data.  The  time- 
averaged  data  (red  line)  shows  a  compression  followed  by  a  sharp  drop  in  pressure  at  a  couple  of 
mesoflap  locations  upstream  of  the  X  =  0  position  and  one  prominent  expansion  followed  by  a 
sharp  compression  near  the  end  of  the  flap  array.  For  the  pattern  noted  upstream  of  X  =  0,  the 
rise  in  pressure  corresponds  to  compression  of  the  fluid  as  it  bends  into  itself  at  the  upstream 
flaps,  which  in  a  time  average  sense  deflect  upward,  followed  by  an  expansion  at  the  trailing 
edge  of  the  flap  to  approach  the  average  pressure  value.  An  expansion  occurs  at  the  farthest 
downstream  flap  as  the  fluid  bends  away  due  to  the  downward  deflection  of  the  flap  (in  the  time- 
average  sense),  and  then  recompresses  to  average  values  at  the  trailing  edge.  The  pressure  data 
does  not  reveal  a  clear  pattern  for  the  remaining  flaps,  which  possibly  indicates  that  these  flaps 
respond  less  dynamically  to  the  pressure  loads. 

The  deflection  history  (20  K  iterations  =~  3  ms)  of  the  flaps  is  shown  in  Figure  54  with 
the  time-average  position  shown  in  red.  The  flaps  are  represented  as  straight  lines  which  can  be 
thought  of  as  a  projection  of  their  upper  surfaces  on  a  plane  normal  to  its  edge,  ft  is  clear  from 
the  figure  that  the  maximum  deflections  (again,  in  a  time-average  sense)  occur  at  the  two  most 
upstream  flaps  and  the  farthest  flap  downstream.  This  corroborates  the  observations  made  above 
for  the  wall  pressure  data.  Figure  55  shows  a  comparison  of  the  time-averaged  deflections  of 
each  mesoflap  as  calculated  from  the  quasi-steady  and  the  dynamic  FS1  models.  Also  shown  is 
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the  interpolated  pressure  loads  (time-  and  span- averaged  for  the  3-D  calculations)  acting  on  the 
flap  surfaces.  The  figure  shows  that,  in  general,  the  mesoflaps  deflect  to  a  greater  extent  in  the 
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Figure  54:  Deflection  history  of  the  six  flaps 
(represented  by  1-D  beams)  with  the  time-averaged 
position  shown  in  red 

dynamic  modeling.  The  net  pressure  load  which  is  based  on  summation  (with  proper  sign)  of 
span-averaged  pressures  on  the  upper  (negative)  and  lower  (positive)  surface  of  the  flaps  shows  a 
similar  trend.  This  is  consistent,  as  higher  loads  would  result  in  greater  amount  of  deflection. 
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Figure  55:  Deflection  and  interpolated  net  pressure  load  (on  flap  surfaces)  comparisons  for  the 
quasi-steady  and  dynamic  modeling  of  the  mesoflaps;  loads  -  symbol,  displacement  -  solid  lines. 
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Also  to  be  noted  is  that  the  loads  acting  on  the  upstream  flaps  are  primarily  positive  (in  the 
upward  direction)  while  that  for  the  downstream  flaps  are  primarily  negative  (in  the  downward 
direction),  ft  is  also  clear  that  magnitude  of  the  net  loads  tend  to  increase  as  one  moves  away 
from  the  center  which  also  corroborates  the  deflection  patterns  observed. 

Figure  56  shows  the  frequency  response  of  the  mesoflap  array  due  to  the  time-varying 
loads  acting  on  it.  This  is  not  a  response  spectrum  in  the  sense  used  in  structural  mechanics  but 
more  of  an  investigation  of  what  frequencies  might  be  playing  a  role  in  determining  the  dynamic 
response  of  the  flaps.  Since  no  damping  has  been  added  to  the  system  and  there  is  no  algorithmic 
damping  provided  by  the  time  integration  method  chosen,  any  initial  transient  or  numerical  noise 
would  not  be  damped  out.  The  reason  for  not  using  damping  was  to  look  into  the  maximum 
possible  response  of  the  system  which  would  have  been  altered  by  the  presence  of  damping.  Also 
since  we  were  not  really  focusing  on  obtaining  a  steady  oscillation  pattern  for  the  mesoflaps,  this 
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Figure  56:  Frequency  analysis  of  flap-response  (blue)  based  on  the  displacement  of  the  edge  of  the 
flap;  red:  natural  frequencies  (four  shown),  black:  low  frequency  shock  motion,  green:  boundary 
layer  turbulence  frequency,  pink:  frequency  based  on  recirculation  time  in  cavity 

was  justified.  As  such  the  peaks  in  the  plot  are  considered  to  represent  the  driving  frequencies 
responsible  for  the  large  scale  flap  response.  There  are  two  predominant  peaks  in  the  plot  for 
each  of  the  flaps  and  both  of  these  are  less  than  the  fundamental  frequency  of  the  mesoflap  (first 
natural  frequency  of  the  modeled  cantilevered  beam).  Some  relevant  frequencies  (and  also  the 
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first  four  natural  frequencies)  are  also  shown  in  the  figure  to  evaluate  whether  any  of  these  have 
an  impact  on  the  flap  dynamics.  The  line  in  black  correspond  to  the  low  frequency  motion 

characteristic  of  similar  shock  /  boundary  layer  interactions,  ~  1 1 5  — ^  [38],  the  green  line 

represents  a  characteristic  frequency  associated  with  the  largest  eddies  in  the  incoming  boundary 
layer,  — ^  and  the  pink  line  is  a  frequency  based  on  a  quantity  termed  here  as  the  recirculation 

time  Trec,  i.e.  the  time  required  for  the  fluid  injected  from  the  cavity  into  the  core  flow  to  re¬ 
enter  the  cavity  and  go  out  again,  and  has  been  derived  as  a  gross  approximation  in  this  work.  To 
derive  this  quantity,  average  near  wall  (parallel  to  wall)  velocities  are  calculated  for  all  the  three 
walls  and  top  (near  flap  array)  of  the  cavity  using  time  and  span-averaged  velocity  data..  A 
further  approximation  is  made  about  the  velocity  at  the  top  wall  (basically  flap  region)  of  the 
cavity  which  is  evaluated  as  the  average  of  the  free-stream  velocity  Ur  and  velocity  calculated 

as  described  above.  rrec  is  then  calculated  as  the  time  it  will  take  for  any  scalar  to  traverse  the 
entire  perimeter  of  the  cavity  based  on  the  average  velocities  computed.  The  frequency  is  then 
given  by  1  /  rrec .  As  evident  from  the  figure,  it  seems  that  the  frequency  which  seems  to  have  the 

strongest  influence  on  the  flap  response  is  the  low  frequency  shock  motion.  The  computational 
time  used  for  the  present  calculations  (after  initial  transients)  is  ~  3.6  ms  which  roughly 

corresponds  to  about  4.45  time  periods  of  the  low  frequency  oscillation.  The  magnitude  of  the 
dynamic  response  matches  well,  qualitatively,  to  the  flap  responses  as  evident  in  Figure  54.  For 
illustration,  the  maximum  amplitude  of  the  response  in  Figure  56  is  for  flap  2  which  can  be  seen 
clearly  by  the  range  of  the  tip  motions  in  Figure  54. 

6.4.6.  Oblique  shock  /  boundary  layer  interaction  with  mesoflap  array  -  data  comparisons 


Figure  57:  Centerline  wall  pressure  for  SBLI  with  control;  pressure  at  the  lower  wall  of  cavity  is 
shown  for  the  region  of  the  mesoflap  array  (-5.6  to  5.6  on  the  X  axis) 
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Figure  57  shows  the  centerline  wall  pressure  data  for  the  experiment  and  computations. 
The  numerical  predictions  are  in  reasonable  accord  with  the  wall  pressure  at  the  upstream  and 
downstream  regions  of  the  cavity  on  the  whole.  The  pressure  rise  starts  further  upstream  in  the 
domain  than  in  the  no-control  case,  which  is  due  to  a  stronger  leading  edge  oblique  shock.  Data 
is  not  available  in  the  region  just  upstream  and  downstream  of  the  cavity  due  to  the  stringer  plate 
(to  which  the  flap  array  is  epoxied).  The  higher  pressure  observed  downstream  of  the  interaction 
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Figure  58:  Streamwise  turbulence  intensity  profiles  at  various  streamwise 
locations;  a)  inflow,  b)  interaction  region,  and  c)  downstream 


is  still  present  in  these  calculations  as  it  was  for  the  no-control  case.  The  pressures  along  the 
lower  wall  of  the  cavity  agree  well  with  the  experimental  data. 

Centerline  axial  velocity  profdes  are  shown  for  two  downstream  stations  in  Figure  58a. 
RANS  predictions  are  in  excellent  agreement  with  the  data,  while  the  LES/RANS  computations 
predict  a  fuller  profile  near  the  wall.  Figure  58b  shows  streamwise  turbulence  intensity  profiles 
at  two  streamwise  locations  downstream  of  the  interaction.  The  computed  turbulence  intensity 
agrees  well  with  experimental  data  at  the  ( A  -  A0 )  /  S0  =8.73  station.  At  the  station  furthest 

downstream,  the  turbulence  intensity  is  over-predicted  in  the  outer  region  of  the  boundary  layer. 
The  peak  turbulence  intensity  at  the  downstream  locations  increases  when  control  is  added,  as 
can  be  seen  by  comparing  Figure  58b  with  Figure  48.  The  combination  of  unsteady  blowing 
and  suction  as  induced  by  the  mesoflap  dynamic  response  appears  to  amplify  turbulence  to  levels 
higher  than  encountered  in  the  case  without  control.  This  effect  is  opposite  to  that  provided  by 
strong,  nominally-steady  bleed  rates,  as  discussed  in  Section  6.3. 
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6.5.  Hyper-X  Boundary  Layer  Trip  Arrays 


A  version  of  the  IB  flow  solver  was  used  by  Meelan  Choudhari  of  NASA  Langley 
Research  Center  to  compute  laminar-flow  base  states  induced  by  some  of  the  boundary-layer  trip 
geometries  used  in  NASA’s  Hyper-X  program.  Trip  geometries  similar  in  shape  to  the  micro 
vortex  generators  were  first  considered  [P-3,  P-4],  followed  by  square-peg  trips. [P-7].  Each  trip 
was  rendered  as  an  immersed  boundary.  The  resulting  converged  base  states  were  then  analyzed 
using  linear  stability  theory  to  determine  the  most  amplified  modes.  Figure  59  shows  the 


Figure  59:  Vortical  flows  induced  by  Hyper-X 
boundary  layer  trips 

formation  of  longitudinal  vortical  structures  in  the  Hyper-X  forebody  boundary  layer  due  to  the 
trip  arrays.  The  reader  is  referred  to  Refs  P-3,  P-4,  and  P-7  for  further  details  of  these 
calculations. 

6.6.  University  of  Michigan  3D  Shock  /  Boundary  Layer  Interaction 

The  immersed-boundary  method  was  also  used  in  a  recent  calculation  focused  on  3D 
shock  /  boundary  layer  experiments  conducted  at  the  University  of  Michigan.  These  experiments 
were  used  in  a  blind  challenge  test  to  ascertain  the  ability  of  current  methods  (mainly  RANS  and 
LES/RANS)  to  predict  the  detailed  features  of  such  interactions.  The  Michigan  team  used  stereo 
PIV  techniques  to  develop  a  reasonably  complete  map  of  the  velocity  field  and  second-moment 
quantities  within  the  interaction  region.  The  shock  generator  was  mounted  in  the  middle  of  their 
wind  tunnel  to  avoid  tunnel-starting  issues.  A  consequence  of  the  mounting  strategy  was  that  the 
shock  generated  was  not  two-dimensional,  but  rather,  conical.  The  entire  interaction  was  also 
influenced  by  sidewall  shock  /  boundary  layer  interactions.  The  NCSU  team  participated  in  the 
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blind  challenge  study  and  was  able  to  complete  one  calculation  of  this  flow  using  the  newest 
version  of  the  LES/RANS  method  [27].  The  immersed-boundary  technique  developed  in  this 
work  was  used  to  model  the  effects  of  the  shock  generator.  The  non  mass-conservative  variant, 
which  reconstructs  the  pressure  field  within  the  band  cells,  was  employed,  as  the  use  of  the 
conservative  form  led  to  band-cell  oscillations  in  pressure  which  then  propagated  along  Mach 
lines  toward  the  surface.  To  account  for  unresolved  boundary  layer  growth  on  the  shock 


Figure  60:  Instantaneous  Mach  number  contours:  University  of 
Michigan  SBLI  experiment 

generator,  the  normal  vectors  defined  at  the  lower  surface  were  transitioned  between  values 
corresponding  to  a  wedge  angle  of  8.5  degrees  to  the  correct  wedge  angle  of  7.75  degrees  over  a 
distance  equal  to  1/10  of  the  length  of  the  shock  generator.  This  modification  was  necessary  to 
obtain  the  correct  shock  impingement  point.  Figure  60  shows  instantaneous  Mach  number 
contours  for  this  calculation.  The  blue  contours  correspond  to  the  blanked-out  portions  within 
the  shock-generator  immersed  body. 

7.  Conclusions 

This  work  has  developed  an  extension  of  an  existing  immersed-boundary  (IB)  method  to 
compressible,  turbulent  flows  and  has  investigated  its  use  in  simulating  the  effects  of  different 
types  of  flow  control  devices:  micro  vortex  generators,  bleed-hole  arrays,  aero-elastically- 
deforming  ‘mesoflaps’,  and  boundary  layer  trips.  The  IB  method  has  been  used  with  RANS 
(Menter  SST)  and  LES/RANS  turbulence  closures.  The  method  has  been  applied  to  oblique 
shock  /  turbulent  boundary  layer  interactions  with  and  without  micro  VG  flow  control 
(Cambridge  University  experiments),  bleed  flow  control  (NASA  Glenn  experiments),  and 
mesoflap  flow  control  (University  of  Illinois  experiments).  The  following  general  conclusions 
may  be  stated  (case-specific  conclusions  may  be  found  in  References  P-1  through  P-9). 

1.  The  LES/RANS  model  originally  developed  in  [15]  provides  good  predictions  of  the 
mean  flow  structure  of  these  interactions  with  and  without  control.  Second-moment 
predictions  are  also  in  good  accord  with  available  experimental  results.  In  general, 
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though,  the  mean-flow  results  are  no  more  accurate  than  those  provided  by  the  Menter 
SST  RANS  model,  which  performs  excellently  for  this  class  of  relatively  weak  shock  / 
boundary  layer  interactions.  Specific  deficiencies  noted  in  the  LES/RANS  model 
include  a  slight  tendency  to  under-predict  the  axial  extent  of  separation  and  a  related 
tendency  to  over-predict  the  recovery  rate  the  boundary  layer  downstream  of  re¬ 
attachment.  The  latter  effect  may  be  overcome  by  a  new  variant  of  the  LES/RANS 
model  [27]  which  is  better  able  to  respond  to  changes  in  the  equilibrium  structure  of  a 
boundary  layer. 

2.  The  compressible,  turbulent  version  of  the  IB  model  performs  very  well  in  mimicking  the 
effects  of  the  control  devices  considered  without  the  need  for  excessive  mesh  refinement 
in  the  vicinity  of  the  devices.  Two  versions  have  been  developed:  one  of  which  is 
rigorously  mass-conserving.  The  non-  mass-conserving  version  provides  smoother  flow 
solutions  near  the  IB.  The  mass-conserving  version  can  result  in  oscillations  in 
thermodynamic  properties  near  the  IB,  but  its  use  is  necessary  for  problems  in  which 
precise  mass  conservation  is  essential.  Examples  of  this  include  flow  through  small  slots 
and  holes. 

3.  The  IB  framework  has  been  extended  to  include  aero-elastic  effects  as  modeled  through 
beam-bending  theory.  Predictions  of  meso-flap  deflection  for  the  University  of  Illinois 
experiment  and  the  associated  effects  on  the  fluid  dynamics  of  the  interaction  appear  to 
be  qualitatively  correct. 

4.  The  Cambridge  University  SBLI  experiments  with  micro  VGs  are  significantly 
influenced  by  three-dimensional  effects  associated  with  sidewall  boundary  layer 
separation.  These  contribute  to  the  enlargement  of  the  primary  zone  of  axial  separation 
near  the  centerline  of  the  wind  tunnel.  A  quantitative  assessment  of  the  effects  of  the 
micro  VGs  is  made  more  difficult  because  of  these  effects,  but  the  general  trend  is  that 
the  vortices  induced  by  the  micro  VGs  deform  the  region  of  separation  but  do  not 
eliminate  it  entirely.  Computations  of  idealized  SBLIs  with  micro  VG  control  show  that 
the  strengths  of  the  generated  vortices  are  maintained  for  a  reasonable  distance 
downstream  of  the  interaction,  and  the  LES/RANS  model  in  particular  predicts  a 
broadening  of  the  vortices  that  energizes  the  entire  recovering  boundary  layer.  These 
trends  suggest  that  the  micro  VG  concept  has  the  potential  of  controlling  weaker  shock  / 
boundary  layer  interactions  (generated  as  part  of  a  shock  train)  even  if  the  primary 
separation  region  is  not  wholly  eliminated. 

5.  The  disturbances  induced  by  bleed-hole  arrays  of  the  type  used  in  the  NASA  Glenn 
experiments  are  highly  localized  to  each  bleed  port  and  are  not  representative  of  a 
continuous  suction  effect.  The  fact  indicates  that  existing  bleed  models  may  have  to 
consider  variable  porosity  effects  for  better  predictions.  Reynolds-stress  predictions 
obtained  from  the  LES/RANS  model  indicate  that  the  use  of  high  levels  of  boundary- 
layer  bleed  suppresses  the  amplification  of  turbulence  often  noted  in  shock  /  boundary 
layer  interactions.  This  result  indicates  that  the  dominant  source  of  turbulence 
amplification  in  un-controlled  interactions  may  be  the  dynamics  of  the  separation  /  re¬ 
attachment  event. 

6.  The  UIUC  mesoflap  concept  as  applied  to  oblique  SBLIs  is  not  a  particularly  effective 
means  of  control  if  the  intent  is  to  remove  regions  of  axial  separation.  Strong  blowing 
into  the  boundary  layer  occurring  at  the  first  few  mesoflap  stations  forces  the  boundary 
layer  away  from  the  surface,  increasing  the  probability  that  the  impinging  shock  will 
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separate  the  boundary  layer.  The  suction  is  applied  near  the  place  where  the  separated 
boundary  layer  would  re-attach  naturally  and  thus  has  little  effect.  Turbulence  levels  in 
the  recovering  boundary  layer  are  enhanced  with  meso-flap  control,  relative  to  the  case 
without  such  control.  This  may  make  the  boundary  layer  more  resistant  to  subsequent 
shock  impingements.  An  analysis  of  the  frequency  response  of  the  mesoflap  system 
indicates  that  the  dominant  forcing  frequency  is  comparable  to  the  low-frequency  modes 
naturally  observed  in  uncontrolled  SBLIs  with  flow  separation.  This  provides  additional 
evidence  that  the  mesoflap  concept  does  not  remove  separation  to  any  significant  extent 
and  may,  in  fact,  increase  separation. 


8.  Directions  for  Future  Work 

Some  challenges  and  directions  for  future  work  are  summarized  as  follows. 

1.  The  IB  method  needs  to  be  improved  so  that  the  best  features  of  the  non-  mass- 
conservative  and  conservative  variants  are  combined.  The  fact  that  the  surface  normal 
vectors  can  be  altered  so  that  they  may  not  exactly  conform  to  the  true  surface  normal  (as 
used  in  Section  6.6)  may  provide  directions  for  future  development,  as  the  normal  vector 
could  be  considered  a  variable  that  may  be  adjusted  to  enable  better  mass  conservation. 

2.  The  coupling  between  the  structural  solver  and  the  IB  movement  algorithm  needs  to  be 
improved.  The  fact  that  mass  is  ‘swept’  from  one  class  of  cell  to  another  as  the  IB  moves 
is  not  directly  accounted  for  in  the  current  implementation.  Techniques  used  to  maintain 
particulate  mass  conservation  in  the  presence  of  IB  body  movement  [39]  were 
implemented  but  were  not  successful. 

3.  Observed  weaknesses  of  the  LES/RANS  method  include  its  tendency  to  under-predict  the 
extent  of  axial  separation,  its  tendency  to  over-predict  the  rate  of  recovery  of  the  re¬ 
attaching  boundary  layer  and  the  lack  of  generality  of  the  calibration  procedure.  A  new 
model  [27]  that  is  designed  to  overcome  some  of  these  deficiencies  is  under  development. 

4.  The  effects  of  three-dimensionality  on  the  flow  structure  of  shock  /  boundary  layer 
interactions  and  the  efficacy  of  various  control  devices  in  these  situations  are  not 
understood  at  present.  A  future  focus  will  be  to  apply  the  LES/RANS/IB  models  to  such 
flows  considering  the  full  effects  of  wind-tunnel  boundary  layers,  sidewall  SBLIs,  and 
comer  vortical  structures.  The  case  presented  in  Section  6.6  is  a  starting  point,  but  much 
more  work,  performed  in  close  collaboration  with  experimental  efforts,  is  needed. 

9.  Publications  Resulting  from  this  Study 

[P-l]  Ghosh,  S.,  Choi,  J-I.,  and  Edwards,  J.R.  “RANS  and  hybrid  LES/RANS  Simulation  of  the 

Effects  of  Micro  Vortex  Generators  using  Immersed  Boundary  Methods”  AIAA  Paper  2008- 

3728,  Presented  at  the  38th  AIAA  Fluid  Dynamics  Conference,  Seattle,  WA,  June,  2008. 

[P-2]  Ghosh,  S.,  Choi,  J-I.,  and  Edwards,  J.R.  “Simulations  of  Shock  /  Boundary  Layer 

Interactions  with  Bleed  using  Immersed  Boundary  Methods,”  AIAA  Paper  2009-1330,  Jan. 

2009.  Presented  at  the  47th  Aerospace  Sciences  Meeting  and  Exhibit,  Orlando,  FL  ,  Jan.  2009. 
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[P-3]  Choudhari,  M.,  Li,  F.,  and  Edwards,  J.R.  “Stability  Analysis  of  Roughness  Array  Wake  in 
High-Speed  Boundary  Layer”  AIAA  Paper  2009-170,  Presented  at  the  47th  Aerospace  Sciences 
Meeting  and  Exhibit,  Orlando,  FL  ,  Jan.  2009 

[P-4]  Choudhari,  M.,  Li,  F.,  Chang,  C.L.,  and  Edwards,  J.R.,  "On  the  Effects  of  Surface 
Roughness  on  Boundary  Layer  Transition,"  Proceedings  of  the  International  Conference  on 
Aerospace  Engineering  and  Exhibition,  2009,  Bangalore,  India,  May  18-22,  2009. 

[P-5]  Ghosh,  S.,  Choi,  J-I.,  and  Edwards,  J.R.  “Numerical  Simulations  of  the  Effects  of  Micro- 
Vortex  Generators  using  Immersed  Boundary  Methods”  AIAA  Journal,  Vol.  48  No.  1,  2010, 
pp. 92-103 

[P-6]  Ghosh,  S.,  Choi,  J.-I.,  and  Edwards,  J.R.  “Simulation  of  Shock  Boundary  Layer 
Interactions  with  Bleed  using  Immersed  Boundary  Methods”  Journal  of  Propulsion  and  Power, 
Vol.  26,  No.  2,  2010  (to  appear) 

[P-7]  Choudhari,  M.,  Li,  F.,  Chang,  C.,  King,  R.,  Edwards,  J.R.,  and  Kegerise,  M.,  “Laminar- 
Turbulent  Transition  Behind  an  Isolated  Roughness  Element  in  a  High-Speed  Boundary  Layer” 
AIAA  Paper  2010-1575,  Presented  at  the  48th  Aerospace  Sciences  Meeting  and  Exhibit,  Orlando, 
FI,  Jan.  2010. 

[P-8]  Edwards,  J.R.  “Hybrid  LES/  RANS  Simulation  of  the  Effects  of  Boundary  Layer  Control 
Devices  Using  Immersed  Boundary  Methods”  Final  Technical  Report,  AFOSR  Grant  FA9550- 
07-1-0191,  Feb.,  2010. 

[P-9]  Ghosh,  S.,  Choi,  J.-I.,  and  Edwards,  J.R.  “Numerical  Simulation  of  the  Effects  of 
Mesoflaps  in  Controlling  Shock  /  Boundary  Layer  Interactions”  AIAA  Paper  2010-4465.  To  be 
presented  at  the  40th  AIAA  Fluid  Dynamics  Conference,  Chicago,  IL,  June,  2010. 

10.  Personnel 

Jack  R.  Edwards,  Principal  Investigator 
Jung-11  Choi,  Research  Assistant  Professor 

Santanu  Ghosh,  Graduate  Research  Assistant  (will  defend  his  Ph.D  Dissertation  entitled 
“Simulations  of  Shock  Boundary  Layer  Interactions  Using  Immersed  Boundary  Methods”  in 
March,  2010) 

11.  Connections  /  Collaborations  /  Invited  Talks 

•  Navier-Stokes  code  with  immersed-boundary  methodology  provided  to  Meelan  Choudhari  of 

NASA  Langley  and  used  to  predict  flows  over  boundary-layer  trip  arrays 

•  Recycling  /  rescaling  turbulence-generation  module  provided  to  Datta  Gaitonde  at  AFRL/RB 

and  implemented  into  AFRL’s  FDL3DI  solver  during  the  Pi’s  visit  to  AFRL  in  July,  2009. 
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•  Invited  talk  on  AFOSR-sponsored  work  entitled  “Hybrid  LES/RANS/Immersed-Boundary 

Methods  for  Simulating  Effects  of  Boundary-Layer  Control  Devices”  presented  at  the  1st 
Shock  Wave  /  Boundary  Layer  Interaction  Workshop  Spring  Meeting,  Ohio  Aerospace 
Institute,  April  15-16,  2008 

•  Invited  talk  on  AFOSR-sponsored  work  entitled  “Simulations  of  Shock  /  Boundary  Layer 

Interactions  with  Bleed  using  Immersed  Boundary  Methods”  presented  at  the  2nd  Shock 
Wave  /  Boundary  Layer  Interaction  Workshop  Spring  Meeting,  Ohio  Aerospace  Institute, 
March  31 -April  1,  2009 

•  Participation  in  blind  challenge  study:  AFRL  Shock  /  Boundary  Layer  Interaction  Workshop, 

Jan  8,  2010 
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